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SUMMARY 


This report demonstrates the successful application of statistical 
variable selection techniques to fit splines. Major, emphasis is given to 
knot selection, but order determination is also discussed. Two FORTRAN 
backward elimination programs using the B^spl ine basis were developed, and 
the one for knot elimination is compared in detail with two other spline- 
fitting methods and several statistical software packages. An example is 
also given for the two-variable case using a tensor product basis, with a 
theoretical discussion of the difficulties of their use. 

1. INTRODUCTION 

Polynomial splines have often been employed in modeling or data fitting 
when the functional form of the relationship between the dependent and inde- 
pendent variables is unknown. The major problem has been how to avoid 
under- or overfitting the data. A strictly mathematical approach is to add 

knots one at a time and move them around until the L 2 (or some other) norm 
of the errors is less than a preselected tolerance level (ref, 1). A major 
problem with this approach is that a good fit depends entirely on the sub- 
jective selection of the tolerance level. A fitting method which attempts 
to avoid this problem is the smoothing technique introduced by Reinsch 
(ref. 2), but it requires the experimenter to have good a prior i information 
about the data or the process which generated it. Both of these methods are 
currently feasible only for functions of a single variable. 

A statistical approach to the curve-fitting problem using the method of 
cross-validation was introduced by Wahba and Wold (ref. 3). The major 



advantage of this procedure is its automation: no a priori information is 

needed. There are several disadvantages, however. Every data point is a 
knot so that the resulting functional form is difficult to use and inter- 
pret, In addition, if there are clearly identifiable trends in certain 
portions of the data such as linearity or sharp bends, this information is 
lost analytically even though it shows up when the spline is plotted. The 
practical use of this technique is also currently restricted to functions of 
one or two variables. The two variable case is considered in Wahba (ref. 

4), with higher dimensions discussed in Wahba and Wendelberger Cref. 5), 

Other statistical approaches to the variable knot spline problem have 
considered the knots as parameters in the model. However, this presents 
problems in finding the least squares solution and in subsequent statistical 
estimation and testing procedures because the model is nonlinear. Tradi- 
tional (ref. 6) as well as Bayesian (ref. 7) approaches have been investi- 
gated, but both are limited in scope and application. Further, in most 
cases, though the knot locat ions have been variable, their number has been 
fixed a priori by the analyst. Some exceptions are the works of Ertel and 
Fowlkes (ref, 8), Smith and Smith (ref. 9), and Agarwal and Studden (ref. 

10), but, as with most other approaches mentioned above, they have not been 
developed to fit splines in several variables. 

The technique investigated in this research is the use of variable 
selection procedures to fit splines. If a pool of knots is fixed in advance, 
then statistical linear models theory can be applied in a variable selec- 
tion framework. There are four major advantages of the variable selection 
approach to fitting splines. First, variable selection procedures are 
essentially user independent (automatic) in their use of the F test as a 
stopping criterion. Second, they are widely available in statistical soft- 
ware. Third, final fits may have straightforward intepretat ions because of 
their simplicity or theoretical foundation. Fourth, regression diagnostics, 
such as out lier detect ion, may be performed. These advantages and other 
desirable properties are discussed in Section 4, along with a comparison of 
several methods and software. 

The theory applies not only Co splines in a single variable, but also 
to splines in several variables using a tensor product basis. However , as 


the careful and detailed developroent of this technique in the one-variable 
case is considered a crucial step to its use in several variables, discus- 
sion of the niult ivar iate case is restricted to Section 7, and includes an 
example of its successful application to aerodynamic modeling. 

The major emphasis of this report is the application of variable selec- 
tion procedures for choosing the number and location of the knots for 
splines in a single variable of fixed order (= degree +1). A detailed dis- 
cussion of this "knot selection" approach is given in Section 2 with exam- 
ples, comparison of methods and software, and applications in Sections 3-5. 
Choosing the spline order with the number and location of the knots fixed is 
of less interest and considered in Section 6 only. FORTRAN programs which 
apply backward elimination in these two contexts’ were written as part of 
this research and discussed in Sections 2 and 6. Their documentation, flow- 
charts, and listings are given in the Appendix. 
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dependent variable 
breakpoints for y 
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aileron, elevator, and rudder deflection 

random error 

mean 

standard deviation 
gridpoint (x^,y J 


knot selection 

mean squared error 

Smith-Smith 

error sum of squares 

Wahba-Wold 


2. THE KNOT SELECTION (KS) PROCEDURE 


Statistical variable selection procedures can be used as a KS procedure 
to choose the number and location of knots in fitting splines. The 
function basis is suitable for this, at least theoretically, because it is. 
easily interpreted. Knots and knot multiplicities correspond to individual 
terms so that selection or deletion of terms is equivalent to selection or 
deletion of knots. The knots are thus selected indirectly. For example, a 
continuous linear spline with knots ^2*’ written as 

£ 

6 +3,x+i P.(x-t.). where u ~u for u > 0 and zero othe r- 

o 1 2 1 1 > 
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wise. Selection of the "spline term" (x - actually selection of the 

knot t.* Because we don't know where the breakpoints should be, we provide 
as candidate variables a liberal number of spline terms, i.e., a pool of 
knots, more than we expect or want to eventually use, and blanket the 
domain. thus, the actual niraber and location of the knots used in the final 
model is unknown at the beginning in the sense that we are selecting from a 
larger set . 


While functions are easily defined in current statistical software 

packages and fit into the statistical hypothesis testing framework without 
modification (ref. 11), computational problems such as carry-over in round- 
off error and mult icol 1 inear ity greatly restrict their use. As will be seen 
in Section 4, the backward elimination (stepdown) procedures are especially 
troublesome because all terms must be fit initially. An alternative is the 
use of the computat ional ly advantageous B-spline basis (ref. 1). Unfortu- 
nately, it does not fit easily into the hypothesis testing framework and 
cannot be used in existing statistical software packages. There was thus a 
need for the development of a KS procedure using B-splines. Construction 
of hypotheses which are useful in B-spline regression, including testing the 
importance of knots, has been detailed in Smith (ref. 12). As part of this 
research, these results have been implemented in two FORTRAN computer 
programs, one of which accommodates the backward elimination of knots using 
the B-spline basis. Examples in Section 3 give the results of using this 
FORTRAN program, and comparisons with several statistical software packages, 
as well as with other statistical spline-fitting methods, are detailed in 
Section 4. 


The use of variable selection is a sort of compromise between the tech- 
niques which use either fixed or variable knots. Its most important advan- 
tage, and one which makes possible all others, is that because the maximum 
number and location of the knots is fixed in advance, the statistical theory 
of general linear models applies. Consequently, the least squares solution 



is easily obtained at any given step, and hypothesis testing and interval 
estiraation are straightforward. As mentioned earlier, details for using the 
B-spline basis are given in reference 12. The selection of knots can thus 
be accomplished through t tests. This fits exactly into the variable 
selection framework for (1) spline models in a single variable, (2) models 
in several variables with spline terms in one or more variables, and (3) 
models in several variables with tensor products defining higher dimensional 
splines. Also, trends in the data in one or more variables may be easily 
detected through the selection of a few knots, Several examples of this 
will be given in the next section. Further, in some experimental situa- 
tions, models may be easily interpreted because the coefficients are physi- 
cally meaningful, as in some examples in Sections 5 and 7. 

3. EXAMPLES OF THE KS PROCEDURE 

Four data sets were examined using the FORTRAN knot elimination pro- 
grOT, The maximum number of continuity constraints allowed for any given 
order were imposed. The first data set, the Indy data, is rather simplistic 
but has appeared in the statistical literature several times in connection 
with curve-fitting with splines. It is a record of the average winning 
speeds at the Indianapolis 500 from 1911-1971, except for 1917-1918 and 
1942-1945, during the two World Wars when the race was not run. Poirier 
(ref, 13) fit the data with a cubic spline with 2 knots, one each at the 
midpoint of the non-racing years. The data were coded so that x - year - 
1910 with knots 7.5 and 33.5. The output and graphs from the knot elimi- 
nation routine are shown in Figures 3.1 to 3.4, with circles around the 
function values of the knots. Using an F-table value of 8.0 (a - 0.01), the 
KS procedure eliminates both knots so that a cubic polynomial is adequate 
to fit the data. If a linear rather than a cubic spline is fit, only the 
knot at x ■ 7.5 can be eliminated (Figures 3.5 to 3.7). 

The second example is noisy data generated from the function used in 
reference 3 

f(x) = 4.26(e - 4e + 3e ) 
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3.1. Output for knot elimination. Indy data. Cubic spline. 
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Figure 3.2. First step of knot elimination. Indy data. Cubic Spline. 



Figure 3.3. Second step of knot elimination. Indy data. Cubic spline 
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Third and final step of knot elimination. Indy data 
spline , 
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for xc[0,3]. For x starting at zero, we generated 100 data points at 
intervals of 1/32 up to 99/32 and added normal random noise, N(vi = 0, 
a * ,2), the value of o the same as that used by Wahba and Wold (WW) . A 
graph of the function and generated data is shown in Figure 3,8, Figures 
3.9 to 3.27 show graphical results of the stepdown procedure for cubic 
splines starting with 19 equally spaced interior knots, and using ah F-table 
value of 8.0. By examining this sequence of graphs, it becomes clear how 
the elimination of knots makes the spline smoother by making it less noise 
dependent. 

An F-table value of 4.0 (a = 0.05) rather than 8,0 results in stepdown 
terminating with 5 knots remaining (Fig, 3.23, p. 20). The latter fit is 
more data dependent and clearly inferior in terms of recovering the desired 
function. Use of the larger F value thus seems appropriate and keeps the 
procedure from terminating **pr emat urely . ” Graphs of starting and ending 
fits to the data, beginning with 39 interior knots, are shown in Figures 
3.28 to 3.29, and the results are roughly the same as when 19 knots are used 
initially (Figure 3.27, p, 22). A phenomenon which occurs throughout most 
of these fits is the downward hook in the upper range of the x*s due to a 
cluster of 3 data points. Figure 3.30 shows the conclusion of stepdown with 
those 3 points omitted and helps to illustrate the fact that different noise 
results in different fits. 

The method used by Wahba and Wold to recover the function is a modifi- 
cation of the smoothing technique introduced by Reinsch (ref. 2). Hiey use 
c ros s^al idat ion to determine the smoothing parameter, and their resulting 
fit is shown in Figure 3,31. Referring again to Figure 3.27, p, 22, we see 
that the results of the two methods compare very favorably. A more detailed 
comparison of these methods and others is made in the next section. 

Smith and Smith (SS) (ref. 9) examine a scaled version of the WW 
function, f(x) = 4.26 [e - 4e + 3e for xe[o,l], A sample 

of size 600 equally spaced points was generated, and a variance of 0.039 (as 
in SS) was used for the normally distributed zero mean noise. Results from 
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Figure 3.9. First step of 
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Figure 3.13. Fifth step of knot elimination. WW data. 
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Figure 3.14. Sixth step of knot elimination. W data. 



Figure 3.15. Seventh step of knot elimination, W data 
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Figure 3.16. Eighth step of knot elimination. W data. 



Figure 3.17. Ninth step of knot elimination, WW data. 
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Figure 3.21. Thirteenth step of knot elimination. WW data. 
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WW DATA 

Figure 3.22. Fourteenth step of knot elimination. WW data. 
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Figure 3.23. Fifteenth step of knot elimination. WW data. 
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two stepdown runs fitting cubic" 
beginning with 19 and 49 knots, 
slightly wigglier fit than that 
These results show that a larger 
possibly a fewer number of final 
simplicity, is more desirable. 


splines are shown in Figures 3.32 to 3.33, 
Three knots remain in Figure 3.32 with a 
in Figure 3.33 with one remaining knot, 
knot selection pool allows reduction to 
knots and a smoother fit, which, for 


Smith and Smith use asymptotic results to determine a stopping rule for 
adding knots one at a time to the model. Figure 3.34 shows their results 
using cubic splines overlaid on the true function. The data were not plot- 
ted so that the distinctions between the two functions would not be lost. 
Applying stepdown using these 9 initial knots resulted in Figure 3,35, a fit 
which smooths the wiggles visible in Figure 3, 34. As seen in the two 
previous figures, however, using a larger pool of knots results in a 
smoother and more satisfactory recovery of the function. The SS method is 
compared in more detail to both the WW and KS methods in the next 
section. 


The final function examined is f(x) = sin (x^) for xe[o,4.5], which 
allows for more than two periods of the sine wave and gradually increases 
the frequency. Three hundred data points were used with a = .2 for the 
normal noise. Beginning and ending cubic spline fits from a stepdown run 
are shown in Figures 3.36 to 3.37, starting with 19 interior knots and 
ending with 9. We note that more knots are needed for the final fit than 
for the functions previously discussed due to the increased curvature of the 
function, >fc>st of the wiggliness in the initial spline fit occurs on the 
more gradual slope at the lower end of the x-range and is removed as knots 
are removed. This phenomenon also occurs on the "flat" portion of the SS 
and WW data. 

In order to assess the effects of a lower noise level on the KS tech- 
n ique , random variables used for the noise on the WW function were generated 
using o = .1 and .05, Final fits are shown in Figures 3.38 to 3.39, and 
referring back to Figure 3.27, p. 22, which shows results using a = .2, we 
see that fitting data with a lower noise level results in more knots 
remaining at the end of the procedure. This tendency is especially striking 
when data from the function itself is fit, that is, when no noise is added 
so that to 
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recover the function we actually need to interpolate. A stepdown from 19 
knots results in a spline with' i2 knots as shown in Figures 3.40 to 3.41, 
Both the true and fitted functions are graphed, but there is no perceptible 
difference between the two, 

Wiggliness in data should be smoothed (i.e., ignored) if it is per- 
ceived as noise, but should be fit if it is perceived as trends in the 
underlying process. Thus, a danger in applying the KS technique is using 
too small or too large a pool of knots. The former problem is illustrated 
quite well in Figures 3.42 to 3.43, where noisy data generated from sin (x^) 
is fit with the KS technique beginning with too few knots to allow the 
bending necessary to recover the function, especially near the third peak. 

It is interesting to see that the three knots eliminated were in the lower 
end of the x range where the underlying function is not wiggly, A 
comparison of Figures 3.37, p. 28, and 3.42 reveals that both have 9 knots, 
but a better fit is obtained from the one which began with 19 knots (Fig. 
3.37): its 9 knots are more selectively and better placed. 

4. COMPARISON OF METHODS AND SOFTWARE 

In the previous section, two functions introduced in the literature (WW 
and SS) were examined using the FORTRAN knot elimination program. The pur- 
pose was to compare results, which we do in this section, in light of what 
we consider to be the most desirable properties of curve-fitting with 
splines. These are; 

( 1) good results ; 

(2) computational efficiency; 

(3) diagnostics capab il i t ies ; 

(4) user independence; 

(5) ease of interpretation; and 

(6) ease of use. 

We also give in this section the results of using several statistical soft- 
ware packages on the Indy and WW data, fitting both linear and cub ic 
splines. 
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Most statisticians have ready access to variable selection procedures, 
either in programs they have written themselves, or in widely available 
statistical software packages. Fitting splines through knot selection with 
these programs is a potential advantage of their use, which is realized only 
if good results are obtained. A summary of the results of using four such 
packages is given in Table 4.1: SAS (ref. 14), SPSS (ref. 15), MINITAB 

(ref. 16), and BMDP (ref. 17). 

Table 4.1. Results of using variable select ion techniques to fit 
splines with four statistical software packages. 


SAS 


Stepwise 


linear 
cub ic 


Indy 

/ 

/ 


WW^ 

/ 

/ 


St epdown 


Indy 

/ 

/ 


WW^ 

/ 

/ 


SPSS 


linear 
cub ic 


/ 

✓ 


/ 

/ 


/ 

/ 


/ 


MINITAB linear 
cubic 


/ 

/ 


/ 

X 


/ 

X 


/ 

X 


BMDP^ linear 
cub ic 


X 

X 


X 

X 


X 

X 


X 

X 


^ Selection pool of 19 interior knots. 

^Numerical output has some inaccuracies, but overall results are correct, 
^Tolerance cannot be made low enough to force entry of necessary terms. 


In the case of stepwise procedures, accuracy was determined by comparing 
outputs for the various packages among themselves, while outputs for the 
stepdown procedures were compared with the FORTRAN B-spline knot elimination 
program. Results are surprisingly good considering the fact that the ”+" 
function basis must be used. Entries marked with an "X” indicate failure to 
produce accurate results or, sometimes, any results at all due to high 
mul t icoll ineari ty in the models or low tolerance, especially in stepdown. 
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The minimum tolerance allowed for SPSS, 10~12^ used to force entry 

of some of the polynoraial terras or to get results in stepdowu^ For BMDP. 
the tolerance of 0.01 for variable selection was not low enough to force 
entry of necessary terms to get results for any of the cases considered. As 
expected, less trouble was had with fewer knots (Indy data), lower degree 
(linear), and simpler models (stepwise). Stepdown gave accurate results in 
several cases even for a large number of knots, but there are limitations. 
For instance, computational problems were encountered by SAS for the cubic 
WW data with 39 knots. The final models determined by stepwise and step- 
down, however, were either identical or very similar. The occasional user 
of splines could thus safely rely on stepwise procedures from one of several 
packages to give good results. 

Table 4.2 compares several spline-fitting methods: Wahba-Wold (WW) , 

Smith-Smith (SS), and knot selection (KS). As the latter method may be 
implemented through several different computer programs, two statistical 
packages and the FORTRAN knot elimination routine are included. All methods 
give good results for the data examined, though as seen in earlier discus- 
sion, care must be taken when using the statistical packages, especially 
for stepdown. Their use of the function makes them computa t ional ly 

inefficient and can cause severe problems. They are handy, however, for the 
occasional user as is the method which is available as an IMSL subroutine 
(ref. 18). The KS techniques depend on setting an a level for the hypo- 
thesis tests and specifying an initial pool of knots but are otherwise user 
independent. The WW method is "completely automatic," while the SS method 
depends on user application of the stopping criterion. The KS approach in 
general produces results which are easier to interpret. 

Results from this section and from Section 3 show that splines 
fit by knot selection recover the underlying functions quite well and 
compare very favorably with the results of Wahba and Wold and improve upon 
those of Smith and Smith. Though somewhat simplistic, the knot selection 
approach provides an alternative to the method of cross-validation and 
offers a great computational savings. In addition, there is the possibility 
of analytic or physical interpretation in many modeling situations, an 
example of which is given in the next section. 
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Table 4.2, Comparison of spl ine~f itt ing techniques 
and software. 


Desirable Properties WW 

Good results / 

Comput at ional efficiency ^ 

Diagnostics capabilities X 

User independence / 

Ease of interpretation X 

Ease of use occasionally / 


^ SAS is available only on IBM-compatible 




Knot Selection 

S£ 

SAS 

SPSS 

B-Spl ines 
FORTRAN 

/ 

/ 


/ 

/ 

X 

X 

/ 

X 

/ 

/ 

/ 

/- 

/- 

/- 


X 

/ 

/ 

/ 

X 

/-I 


X 


machines . 


5, SOIE SPECIAL APPLICATIONS 

Probably the most useful application of the ^ technique is data- 
sraoothing, and in Section 3 we saw several examples of recovering underlying 
functions from noisy data. A variation that is useful in simulation 
experiments is smoothing the sample quantile function. This 
function is a left-continuous step function defined as Q(u) = 

(i-l)/n < u < i/n, where n is the sample size and i-th 

order statistic. Experimental conditions can be simulated by generating 
data which behaves like the original, and a smoothed sample quantile 
function provides a continuous distribution from which to draw the simulated 
data. An advantage of smoothing the sample quantile function, rather than 
its pseudo-inverse, the sample cumulative distribution function, is that the 
former always has domain [o,l] regardless of the type of distribution. 
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Programming can thus be standardized, as, for example, in the determination 
of the original knot selection pool. 

The KS technique is also useful in modeling. For example, stepwise 
regression has been applied successfully by Klein, Batterson, and Smith 
(ref, 19) to model flight data using splines. They use function terms 

defined in the angle-of-attack variable in a Taylor series expansion of 
force and moment coefficients in order to model longitudinal motion of an 
airplane. One of their simple *'spl ine-mod if ied" Taylor series expansions of 
the vertical aerodynamic force coefficient is given by 


C 

z 


= C (a) , C (a) q* + C (a) 
e 


z q * =0 
6 =0 


6 

e 


where 


C (a) » C (a = 0) + C a + S A. (a - a„ ) 
z z z a ^ ^ ^ ^ 

a i-2 


C (a) = C + 
z z 

q q 


U2 

Z 

£=2 


Bj^(a - a^) 


0 

+ 






0 

+ 


e e 

and a is the angle of attack, q* is the nondiraens iona 1 p^ch rate, 6^ 

is the elevator deflection, C = /9a, C = 3c /^q*, C *= 9C /9 5 , 

z z z z Zx» ze 

a q 6 

e 

They then use stepwise :regression to select terms, and thus knots, in the 
model. This spline representation preserves the concept of stability and 
control derivatives inherent in the usual Taylor series expansion of aero-- 
dynamic coefficients but has the advantage of providing a representation of 
over an extended range of the angle of attack a. A global model over 
the observed range of a is thus obtained through the use of splines. 


Ill 1 
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6, OTHER USEy OF VARIABLE SELECTION PROCEDURES IN SPLINE REGRESSION 


Thus far we have emphasized the use of variable selection to choose the 
number and location of knots. There are other possible, but perhaps less . 
useful, "extensions" to spline regression of variable selection procedures 
based on polynomial or multiple regression models. In the latter cases, the 
purpose is to determine the polynomial degree and the important independent 
variables and interactions. This is accomplished by examining the 
contribution of individual terms in the. model. With univariate spline 
models, however, there are several polynomial pieces, not just one, whose 
degrees may be examined, and, as seen previously, we may examine the 
importance of each knot. Also, one may wish to examine the continuity 
conditions at one or more breakpoints as in the example discussed by Smith 
(ref. 11), Thus, the complexity of the spline model over the polynomial 
model manifests itself in the greater number of ways the dimension of the 
spline parameter space may be altered. Splines in several variables present 
even more possible diversity since, for example, two*-var iable spline 
continuity occurs not across points but along lines connecting grid 
po ints. 

While it might be nice to have a single software package which could 
perform any combination of these spline hypothesis tests, it is neither 
feasible nor desirable. The major reason is that variable order splines, 
i.e., splines with polynomial pieces of different degrees, have not been 
suff ic iently researched by mathematicians to allow for the satisfactory 
construction in a general framework of a basis using either functions or 

S“splines. Lowering or raising the degree of a single polynomial piece must 
be accomplished by applying restrictions to the model, and hypothesis tests 
must then use restricted least squares. In simple cases this may be 
straightforward (references 11 and 20), but in general the task is unmanage- 
able. For example, the user is subject to hidden analytical errors as when 
the regression or hypothesis degrees of freedom are not equal to the number 
of restrictions because some restrictions are obtained automatically through 
linear combinations of others. While theoretically such dependencies can be 
checked., the usual methods would need some revision in the case 
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of B-spline regression since hypotheses involve values of the fitted spline 
or its derivatives (ref. 12). In the case of the function basis, most, 

but not all, of Che individual terms are meaningful. However, the innocent 
yet indiscr iminant selection or removal of terms through hypothesis tests 
can result in fits which are statistically valid yet nonsensical because 
they are uninter prefab le in terms of polynomial degree or knot locations 
(ref. 11). Because of these various difficulties, it is reasonable to con- 
struct task-specific procedures. 

The application of variable selection to knot selection, as in the 
examples in Section 3, is useful for smoothing data with a fixed order 
spline with maximum continuity conditions. In these cases the interest is 
not in the spline order but rather in determining the minimal number of 
knots deemed adequate to faithfully represent the data. Cubic splines are 
popular because of their low degree and second derivative continuity. The 
selective use of forward or backward algorithms in some statistical software 
packages using "+^' functions (see Section 4), or the backward elimination 
FORTRAN program developed here using B-splines, may be used for this 
purpose. 


Another possible "extension” of variable selection to splines is the 
determination of the polynomial degree while keeping the number and location 
of knots fixed, chat is, not consider the knots as "variables" to be either 
entered or removed. Because of the difficulties with variable order splines 
discussed above, we must restrict ourselves Co polynomial pieces of the same 
degree. Unfortunately, even further constraints are necessary for this 
version. The ideal situation would be to compare a maximally continuous 

r u . « ■ . . 

(C ^ k-th order spline, i.e., a k-th order spline with continuous f, 

(1) fk-2) . ^ ^ k-3, 

^ ,...r , with a maximally continuous k-l-st order spline (C ]. 

A formal test, however, is not possible. This can be easily seen by consid- 
ering a specific example using the partial ordering of some spline models 
given in reference 11, Basis elements for and quadrat ic , s plines 

and for linear splines with one knot are shown in Fig. 6.1. A compari- 

son of orders 3 and 2 (degrees 2 and 1) which retained maximum continuity 
conditions would require comparing the quadratic with the linear. 
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Figure 6.1. A partial ordering of some spline spaces. 

Neither is a subspace of the other, however, so they cannot be formally 
compared (via testing). A solution of a sort is available if the 
quadratic and the cP linear are compared, since, as can be seen from the 
figure, the c9 linear basis generates a subspace of the quadratic 

space. 

In general, a test to compare spline orders can be made between splines 

of order k and k-1, both having continuity C . In the case ot cubic 
splines, for example, we could allow continuity of the function and its 
first (but not second) derivative in order to determine whether the order 
could be reduced from 4 to 3 or increased from 3 to 4, Since a C cubic 
has sufficient smoothness (at least to the eye), the procedure is not so 
objectionable. Considerably less satisfactory, however, are the cases for 
linear and quadrat ic spl ines . In comparing splines of order 3 and 2 as. seen 
in Figure 6.1, the quadratic spline would be continuous but not its first 
derivative while in comparing splines of order 1 and 2, the linear spline 
would not even be continuous. . Of course, the results of formal tests can be 
used in combination with informal comparison between SSE*s of the models of 
interest to decide upon an acceptable model, and we recommend this approach, 

A backward elimination FORTRAN program using B-splines has been devel- 
oped for the purpose of reducing spline order using the nesting of some 
'*sub-opt iraal” spaces as described above. Details for the appropriate B- 
spline hypothesis tests are given in reference 12. The listing, documen- 
tation and flowchart for the program are given in the Appendix, and we 
illustrate its use with the Indy data. While some statistical software 
packages could undoubtedly be used by defining functions as in knot 
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SBlcctioHj no 3tC6ropt was rnadc to us6 th6in in this cont6xt. However, usin^- 
the results of Section 4 as a guide, we surmise chat several backward elimi- 
nation procedures wuld be suspect while most forward selection algorithms 
should give fairly accurate results. Again, tolerance levels may have to be 
made small in order to force entry of certain terms. 

Figure 6.2 gives the FORTRAN program output for stepdown order selec- 
tion on the Indy data starting with a cubic spline (order 4) with the two 
knots in rnid-WWI and WII as in Section 3. The program compares splines of 
different order with the same continuity conditions, though other fits are 
given for information purposes. For this example, order reduction is made 
f r om cub ic to Quadratic to. 1 in ear, Es, t ima tes of the B—spline coefficients 
and their standard errors are given for the spline of lowest order which can 
adequately fit the data, and the highest continuity conditions are imposed. 
For this case it is the linear, 

A graphical display of these results is quite helpful, and Figure 6,3 
shows a partial ordering of the relevant spline spaces along with hypothesis 
test results and SSE^s from the program. The dotted lines indicate the 
stepdown caparisons we wish to make, while the solid lines indicate those 
we can actual ly make through formal comparisons (tests). The importance of 
user input into the variable selection process is becoming more widely 
recognized, and here especially, because the formal tests available are not 
exactly what we would like. Consequently, we recomraend the use not only of 
the formal tests, but also of informal comparisons between SSE’s (or MSB's) 
of competing models using a display such as Figure 6,3. 

We illustrate this technique by going through Figure 6.3 step by step, 
and we will discover some interesting characteristics of splines along the 
way. We first observe that while a formal test is not possible between the 
C? cubic and the quadratic, it would not even be necessary since the 

C* quadratic has a smaller SSE than the cubic. A better fit is thus 

obtained with a lower degree! This phenomenon could never happen with poly- 
notnials, but such are the vagaries of splines. An informal comparison in 
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TVC 9100TTCST SPLihC CT OPDCR K- 4 WITH MPlXlrUl COriTiNUlTY C 2 
HfiS SSE* 385.29e0821B fiND ^GE- 7.08306290 

CPh ORDER K- 4 WITV! SUBHiftXirt /1 CONTIMilTY C I 

BE REDUCED TO ORDER K- 3 WITH COfrriMJTlY C 1 7 


YES. 

FOR K* 4 PHD C 1 

FTABLE VPLUE • 4.00000000 OBSERVED F- 1.27322533 

5SE- 340,34966200 rSE- 7,4:169494 


TVE ShOOTfCST SR-l^€ OT ORDER K- 3 WITH nAXlrtil CONTINUITY C 1 
HftS SSE- 376.65994630 PND MSE- 7.53319093 

CPN ORDER K- 3 WITH SUB-mXirtXI CONTIhiJITY C 0 

BE REDUCED TO ORDER K- 2 WITH rtPXiMlI CONTINUTIY C 0 


YES. 

FOR K- 3 PND C 0 

FTABLE VPLUE ■ 4.00000000 OBSERVED F- 3.69172563 

SSE- 366.45371260 HSE- 7.67611901 


* * ¥ ♦♦3rMriM[WWc*2WoMO|nMo|oMcWO*^Mc*3toUu«3r^ 

T>C SNOOTVCST ^IhE Cf ORDER K- 2 WITH flf^.IMLti CONTINUITY C 0 
HPS SSE- 453.46808846 PND fiSE- 8.09153115 

CPN ORDER k- 2 WITH SUB-MAXlrUM CONTINUITY C-1 

BE REDUCED TO ORDER K- 1 WITH MPXlMJi CONTIMJTIY C-1 ’ 


NO. 

FOR K- 2 PND C-l 

FTPBLE VPLUE - 4.00OC00O0 OBSERVED F= 285.46925596 

SSE* 328.52409313 flSE- 6.70457333 


PROCEDURE TERMINPTE5. 


PROCEDURE TERMINPTES 

N coer 

1 73.61675803 

2 88.49577629 

3 114.97961700 

4 157.47B461B2 


WITH L- 3; K- 2; C 0 

ST, ERR. 
2,19886373 
1.11724005 
.94655257 
1.10901397 


FURTVCR IhTORMPTION 
FOR K* 1 PND C-1 

SSE- 6070.37277252 fiSE- 116.73793793 


Figure 6.2, Output for order reductiion. 


Indy data. 
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C ^cub ic 
348 



cubic 


385 


quadratic linear 


quadratic 

368 


C“ ^ linear 
328 


/ 


/ 


not 

'sig. 



•sig. 


C ^ constant 


376 


453 


6070 


Figure 6.3. Partial ordering of spline spaces including SSE's 
and results of order reduction tests. Indy data. 

in going from the quadratic to the linear reveals an increase of 

77 in the SSE. While this increase cannot be formally judged insignificant, 
we may wish to draw such a conclusion based on the results of the formal 
test which compares the quadratic with the C° linear: the larger 

increase of 85 is insignificant in that case. Having thus "safely" arrived 
at the linear, we must decide whether to further lower the order. The 

very large F value (285) from the program output which compares the C“ ‘ 
linear and the C~ ^ constant splines reveals the importance of the linear 
trend. The big increase of 5742 in SSE from the linear to the C~ ^ 

constant is thus highly significant, and since the increase of 5617 from the 
linear to the C~ ^ constant (the desired comparison) is only slightly 
smaller, we conclude that the use of a constant spline fit is inadvisable. 


7. SPLINES IN SEVERAL VARIABLES 


A mathematical theory for splines in several variables is still devel- 
oping, and a "satisfactory" basis even in two variables has not been found. 
However, tensor products of either "+" functions or B-splines can be used to 
form a spline basis in several variables. While a tensor product basis is 
somewhat clumsy and its interpretation difficult, we explain here some theo- 
retical aspects of its use for the two variable case and give an example. 
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As in the example in Section 5', a spline-modified Taylor series expan- 
sion can be used to model aerodynamic force and moment coefficients, TViis 
time, however, Klein and Batterson (ref. 21) use splines in two variables, 
the angle of attack a and the sideslip angle b, to approximate the lat- 
eral force coefficient and the rolling and yawing moment coefficients. They 
use the yawing moment coefficient as a typical example, and Cn 

can be expressed as 

, ^ + C (a)p’ + C r' 

0 * 0 n n 

r p r 

r' = C 

(7.1) 

+ C fctl 6 + C (a) 5 

ng ^ a 

a r 


C = C (a.b). 
n n 0 = 

a 

■ P' = 


where p and r are Che rolling and yawing velocity and and Sj- 

are the aileron and rudder deflection. They approximate the function Cp^(a,b) 

by 


C (a,b) 
n 


C + C,b + 
o 1 


£1 

Z 

i=l 


(A . + A, . b) (a - a . 

01 li 1 + 


+ 



B . (b 
oj 



£l £2 
E E 


i=l j=l 


D . .(b - b .) 
ij J + 


(a 


a.)° 


(7.2) 


while the remaining functions in (7,1) are approximated by splines in a 
alone. Results from a stepwise regression using these terms are not as good 
as in the one-variable case, and some fine-tuning remains. 

From a theoretical point of view, the tensor product basis does not 
have the nice interpretation of knots and continuity constraints as in the 
one-variable case, even using functions. There is, however, a one-to- 

one correspondence between two-variable *‘*<*** function terms and grid points, 
or nodes, and for this reason, we use the term node . basis to refer to tensor 
products of the function basis. As before, we use right-continuous 

functions so that 0^ is 1. Tensor products of B-splines may also be used to 
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construct a basis for splines in two variables, and we shall see that the 
same advantages and disadvantages of the one-variable case carry over. 

We discuss the simplest two-variable case in some detail: first order 

splines, i.e. step functions. Their applicat ion is somewhat limited, but 
there are several reasons for their detailed consideration. First and fore- 
most, splines in two variables are difficult to envision and manipulate, and 

consideration of the simplest case, namely constants, is thus highly desira- 

ble. Second, as seen in the example above and in Section 5, the estimation 
of aerodynamic force and moment coefficients using a spline-modified Taylor 
series expansion reveals the importance of using constants from both inter- 
pretative and numerical points of view. Finally, the two-dimensional cumu- 
lative distribution function is a first order spline in two variables. 

Thus, the constant case i while limited, has already shown its usefulness. 

We first discuss the node basis by way of example. Suppose breakpoints 
in the x variable occur at x^, X 2 , X 3 and in the y variable at y^ 
and y 2 for data in Xq < x < and y^ < y < y^. A function basis 

of order 1 in the x variable is (x - (x - X 3 )^ and in the y 

variable is (y - (y - y 2 ^^-* tensor product basis is formed 

by taking all the 4 x 3 = 12 products (x - (y ” ^ j ^ i = 0 ,..., 3 ; 

j - 0,..., 2. Each basis element in the two variables is thus a plane of 
height one bounded below by the line y = yj and on the left by the line x 
= Its support is thus a quadrant of a sort (lA) . We call the 

intersection of these boundary lines, the corner of the quadrant, a node, 
denoted *ij. Figure 7.1 shows the relevant grid and nodes. Through any 

yz 
yi 
yo 

Xj, Xj X 2 Xg Xi, 

Figure 7.1. Nodes for a tensor product of functions. 
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variable selection procedure, a model may be found whose terras are a subset 
of the 12 basis elements. Such a selection might result, for example, in 
the nodes shown in Figure 7.2 with the statistical model 


f(x,y) = 6 ocroyo+ + ^02X0+72^. + 611 X 471 ^ 


+ 621X2^71+ + 022X2+72 


where x ^ . 

L + -" j + 

the application 


is an abbreviation 
of this technique 


for 

to 


+ + 632x3+72+ + e, 

(x - X.) ° (7 - 7j)+- 
aerodynamic modeling. 


We saw earlier 



Figure 7,2. Nodes resulting after variable selection on 
a tensor product of functions. 


For splines of higher order, Che same principles apply in forming the 
basis elements: they are the tensor product of one-variable functions. 

Knot multiplicities in one variable result in node multiplicities in several 
variables. The absence or presence of a node or node multiplicity corres- 
ponds to the absence or presence of a certain basis element. There is thus 
some carry-over frora the one-variable case in interpreting the role Chat 
basis elements play, and also in the fact that standard variable selection 
software may be used. The major drawback of this basis, as in the one- 
variable case, is computational. The basis elements do not have small 
support, so that roundoff errors get worse as computations increase. 


The computational difficulties present in the node basis lead to con- 
sideration of tensor product B-splines. While Che formulation of the latter 
basis is straightforward, its interpretation and use in model selection 
through hypothesis tests are not. The polynomial degree and importance of 
knots in mode ling are considerations that carry over from one to several 
variables, and unfortunately, so do their difficulties when using B-splines, 
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To compare differences in the two-variable case between the node basis 
and B-spline basis, we consider a simple grid with nodes indicated (*) in 
Figure 7.3, 


— 

2 

c • 

4 

1 

3 


xo xj X2 

Figure 7.3. Nodes for model (7.3), 

The statistical model for first order splines is thus 

f(x,y) = 6 ooxo+yo+ + 001X0^.71^ + + e. (7.3) 


The function is a "true” spline in both variables except when Yi^j 

for then f is constant over fxQ,X 2 ). If this model is represented with 
B-splines, each cell i is the support of a right-continuous plane which 
has height 1. Using Che notation B^(x,y) for Che basis element for each 
cell the model may be written 

4 

f(x,y) = [ 6 B.(x,y) + e subiect to Bi = 63 . 

i = l 

This B~spline model is somewhat more complicated than the function basis 

in its representation because of the model restrictions. it is also not 
obvious how to interpret the B-spline coefficients in terms of the presence 
or absence of nodes. 

These simple examples illustrate that the function terms are iden- 

tifiable and meaningful on a grid as nodes, just as they correspond to knots 
in the one— variable case. They thus hold an advantage over the tensor 



product of B-splines from an interpretative point of view. As B-splines 
hold the computational edge, however, it would be desirable to identify the 
linear combinations of B-splines which correspond to the presence or absence 
of nodes. The interpretation and use of tensor product splines of higher 
order is more difficult and remains to be examined in detail. 

8. SUGGESTIONS FOR FURTHER WORK 

There are potential research areas for both the un iv aria te and multi- 
variate cases. In the univariate case, an efficient stepwise computer rou- 
tine using B"splines could be developed. This would give Che user the 
choice of forward and backward procedures with a computationally efficient 
basis. The use of knot selection to fit data with loops could be investi- 
gated, and approaching the problem using the parametric technique of Smith, 
Price, and Howser (ref. 22), seems feasible. The successful use of splines 
in two variables has already been demonstrated (Section 7), but further work 
remains such as investigating fits to known underlying functions like we 
have done in the one-^var iable case. IVo-d imensional pictures in this case 
would be most helpful. Also, while the multivariate mathematical theory is 
still developing, interpretation of Censor-product bases from a statistical 
perspective could continue from that begun in Section 7. 
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APPENDIX 

Program Documentation 

Two FORTRAN programs have been written which adapt stepdown procedures 
to B-spline regression. One program is for knot elimination while the other 
is for reducing the spline order. Theoretical details and appropriate ref- 
erences are given in Sections 2 and 6. The . programs are written in FORTRAN 
5 and have been implemented on both the ODU DEC-IO and the NASA/Langley CDC 
Cyber computers. Notation is patterned after that of de Boor (ref. 1), and 
definitions of parameters are given in the subroutine VL2NT, the second 
subroutine called. All necessary input is read in or specified in subrou- 
tine DATl: the data, sample size NDATA, (initial) spline order K = degree 
+1, (initial) interior breakpoints and endpoints BREAK(*), number of conti- 
nuity conditions V(») at the breakpoints, number of intervals L = interior 
breakpoints +1 > and tabled F value to be used in hypothesis tests. For 
equal spacing, the breakpoints and continuity conditions are most easily 
specified through a DO loop. Variables are dimensioned by one of three 
parameters (defined in comment statements) which are specified in the PARA- 
METER statement at the beginning of the main program. 

Data must be interior to [BREAK(l), BREAK(L-Hl) ] • For the Indy data, 

"X min = 1 and X max = 61, so we arbitrarily set BREAK(l) = 0 and BREAK(L+1) 

=» 62. V(I) is the number of continuity constraints at BREAK(l). For 
example, V(l) = 0 means that the spline is discontinous at BREAK(l) while 
V(2) = 3 means there are 3 contiguous continuity conditions on the spline f 
at BREAK(2), i.e., f, f', and f” are all continuous at BREAK(2). Note that 
V(I) must be less than or equal to K-1 in order to have a "true” spline, 
not a polynomial, across BREAK(I). We always set V(l) = 0, though only for 

"symmetry" in the endpoint conditions, and V(L+1) need not be specified 

since it is never used nor referred to. 

The subroutine FLAG is designed to catch user input errors which would 
otherwise cause the program to terminate abnormally or give inaccurate 
results which may or may not be obvious to the user. Sample output detect- 
ing errors in the input information of the Indy data is shown in Figure 

A.l. 
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TT€ OraiEP K • 4 


UC ^ INTERVfiLS L - 3 


7>c DiroeiON N • 

5 



BREAKPOINTS 

COKriMJlTV CONDITIONS 

S.00000000 


0 


33.50000000 


4 


7.50060000 
6£. 00000000 


3 


T 

INDE^ 



5.00000000 

1 



5.00000000 

2 



5.00000000 

3 



5.00000000 

4 

KEND( 

D- 4 

7.50000000 

5 

,KEND( 

2)" 4 

6Z. 00000000 

6 



62.09009000 

7 



62.00000000 

9 



62.00000000 

9 




BREAKPOINTS MST BE STRICTLY INCREPSING. 

BREAKPOINT 33.50000000 IS NOT (.CSS THAN BREAKPOINT 

TVC NLMBER OT CONTINUITY CONDITIONS MUST BE STRICTLY 
LESS THAN TVE SPLirC ORDER K. V( 2)- 4 
AT BREAKPOINT 33.50000060 IS TOO LARGE. 

X VALUE OUT OF RANGE. 

X( U> 1.00000000 IS NOT IN T>C RANGE BREAK(1)« 

TO BRCAK(LAST)- 62.00000000 


STEPDOWS CAT^T proceed. PROGRAM ABORTS. 


Figure A.l. Sainple output detecting input errors. 


7.50000000 


5.00000000 


Indy data. 
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Several lines in the programs are for plotting only. These are calls 
to the CDC system subroutines PSEUDO, INFPPLT, and CALPLT and the DO loop 10 
which calculates the spline values at the knots. 

For the knot elimination routine, input data and subsequently 
calculated information are printed by means of subroutines DATl and OUTNTS. 
This includes data values, spline order, number of intervals, dimension of 
the spline space, and knots. At each step of the procedure as indicated by 
the number of intervals L, the F-ratios for the importance of each 
breakpoint are given along with the SSE and MSE. If a breakpoint can be 
eliminated, it is specified and the procedure continues to stepdown. If no 
breakpoint can be eliminated, the resulting number of intervals and spline 
order are given along with a list of the values of the B-spline coefficients 
and their standard errors. Sample output appears in Figure 3.1, p. 8 , in 
Section 3, 

As in the knot elimination program, the subroutines DATl and OUTNTS of 

the order reduction routine print input data and subsequently calculated 

information. In addition, at each step, the printout gives the SSE's and 

K“2 

MSE's for two splines of order K, one with continuity C and the other 

K— 3 

with continuity C . The hypothesis test is described in words with the 
results of the F test indicated. When further order reduction is not 
possible, estimates of the B-spline coefficients and their standard errors 
are given for the spline of lowest acceptable order with highest continuity 
imposed. Additional information is given by including the SSE and MSE of 
the next lowest order spline. Sample output for the Indy data appears in 
Figure 6.2, p . 41 , 

Flowcharts are given in Figures A. 2 to A. 3 followed by the program 
listings. A full listing of the knot elimination program from a CDC Cyber is 
given, including the subrout ines , of de Boor (ref. 1) that are used. For the 
order reduction program we list only the main program and the subroutine 
SSHYP2, a variation of SSHYP appearing in the first program. 



Figure A, 2. Flowchart for knot elimination program. 
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oc fS 


ORiGfWAL PA 

■ i quality 


*-’Ac:g ;:? 


00^. 


'^S'i 


KNOT ELIMINATION PROGRAM LISTING 


PROGWM >«.OT(IH>UT,OUl1^,TflPE6-OUTPin-,TflPE2B.TflP€21.T«»E22) 
C STEPDO^ FOR BRErtO*OINT ELIMINATION FOR FI>€D ORDER K, 

C T>C FUNCTION fK> ITS FIRST K-2 DERIVATIVES MUST BE CONTIHXXS. 

C 

C NDNRK IS AT LEAST TVC SATPLE SI2E, rCATA. 

C WAX IS AT least N. WITH MAXIMJM CONTIMIITY CONDITIONS, 

C N-L-H<-1. UITW NO CONTIHJITY CONDITIONS, N-L»«. 

C KTWAX IS AT LEAST RjRi. 

C 

C 

PARAMETERC rm<- 100, NDMAft<-200, K7rm<-2000) 

REAL BCOEFCWAX), QCKTWAX), DIAGOOWAX), T(NDMAe<) 

* ,DCOEF(rm<), BRTCWAX), BLFCWAX),F(NDMAX) 

» ,CTRAST(rm<),AA(WAX,WAX), ERROR(NDMAN) 

» , MSE, rfiH, SE<WAX),FRATIO<WAX) 

• ,BB(WAX,MtW), LlNV(WAX,N1A>0,FB(Wf«) 

INTEGER ERRDF,HDr,V,KE>0(rf1AX) 

COrtlON /DATA/ NDATA, X(NDMAX), Y(NDNAX}, FTABLE 
COWON /APPROX/ BREAK(NIAX), C0EF(KTN«O, L, K, VCWAX) 

ICOLNT-0 
C ENTER DATA 

CALL DATKICOLHT) 

C GET TVe KNOT SEQUENCE 

CALL VL2NT(BREAK,L.K,V,T,N,KEND) 

C PRELIMINARY OUTPUT 

CALL OUTNTS<BREPK,V,L,T,N,K,KETfl» 

OCCK IhRUT DATA 
IFLAC-0 

CALL FLAG(1FLAG,N) 

IFdFUtt .EQ. 1) GO TO 25 

CPU- PSEUDO 

C TEST FOR CONTIHJOUS K-l-ST DERIVATIVE AT EACH KNOT 
JDERIV-K-l 
I FMIH-FTABLE 

LMl-L-1 

C GET TVC LEAST SQUARES FIT, I.E., TVC B-SPLlrC COEFTICIENTS 
CALL LSTS01(T,N,K,0,DIAG,BCC£n 
C LSTSQl CALLS BSPLV8, BChFfiC, AND BCHSLV 

C GET SSE AND MSE 
ERRDF-NDATA-N 

CALL BSPLPP(T,BCOeF,N,K,DIAG,EREAk,COEF,L) 

CALL ERRL2KF, ERROR. ERRDF, SSE, MSE) 
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C ERR^ CflLLS PPVfiLU 1>HICH CflLLS IN7ERV 

LPl-L+l . 

DO 10 I-l.LPl 

FB( I) -PPVfiOJC BREPK, C0EF,L, K, BREPK( I ) , 0) 
10 C0NTI^^JE 


PLOTS 

CA-L II^OPtT(0,NDATA,X,l,F,l,0,,62,,74.,lS8.,l.,9, 

« 9HINDY DPTfl,l,lHY,0,5,,4,,,75..75) 

CPLL irCC«.TC0,NDPTP,X, 1, Y, 1,0. ,62. ,74. , 15Q. , 1 . ,9, 

♦ 9HINDY DPTP,l,lKf,22,5.,4.,.75,.75) 

CPU- irrOPLTa,U11,BREPK(2),l,rB(2),l,0.,62.,74.,150.,l.,9. 

# SHINDY DPTP,1,1HY,1,5.,4.,,75,.75) 


irCL .^€. 1) GO TO 12 
I4?ITE(20,11) SSE,ttSE 

11 rORMPK//' SSE-',ri6.8,5X, 'NSE-',ri6.0) 


GO TO 9 

C TEST irfWTPMCE OF EACH BREPKPOINT 
12 i^ITE(20,2) UFTPBL£,SSE,rSE 

2 rORflPT(///' L-M3,5X, T-TPBLE VPLUE IS^F16.8,//' 

* * SSE*',F16.8,5X, 'MSE-^F16.B// 

m ' F-ftPTIOS PRE: BREPkPOINTS PRE') 

3 DO 5 II-2,L 

ID-II 

CPU- CNTPSTdD, JDERIV,N,k:,L,T,BR£PK,KEND,BRT,BLF,DCOEr 
» ,CTRPST) 

C CNTRST CPLLS BCOHT 

CPLL S5ftT3(BC0GF,CTRJ=eT, Q,K,H,BRT, VPR,SSH, r«H,HDF) 

C SSHVP CPLLS FORSUB 

FRPTIOan-NSH^MSE 
H?ITEC20,4) FRPTIO(II),BREPK(II) 

4 FORWT<2F16.0) 

IF(FRPTIOai) .GE. FMIN) GO TO 5 
FNIN-FRPTIOdl) 

KNOT* 1 1 

5 CONTINUE 

IF (FMIN -LT. "rrPBLE) GO .TO 7 
i4?rTE(20,6) 

6 F0RhPT</' NO BREPKPOINT CPN BE ELIMINPTEDM 

GO TO 9 


7 FMIN-rRPTIO(KNOT) 
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UPITE<20,0) BRCPKCKNOT) 

8 rORMP^K//' BRC»<P0INT\r7.3, ^ IS Q-IMINPTU)^) 


C RELABEL KNOT SEQUENCE T AS UCLL AS BREAK, KEND, AND V 
CALL REKNOT(KEND,KNOT,N,K,L,T,V, BREAK) 

GO TO 1 

C PRINT RESULTING COETICIENTS, STANDARD CEVIATIONS, AND T-VALUES 
9 CAUL STI€RR(O,0COEr,K,N,L,rtSE,DIAG,AA,SE,LIW) 

C STDEPR CALLS BCHINV AND MAPvEC, 


CALL CALPLT(0.,0.,999) 


25 STOP 
END 

C INDV DATA 

suBROurrirc datkicount) 

COftlON STATDCNTS /DATA/ AND /APPROX/ ARE USED, 

C 

C TNIS SUBROUTirC READS IN THE DATA AND GIVES THE NLHBER AND 
C PLACETENT OF T>€ KNOTS FOR TV€ FITTED SPLirC. 

PARAMETER (rflAX-100, NDMAX-200, K7M1AX-2000) 

INTEGER V 
REAL Y X 

COmON / DATA / NDATA, X(NDhAX), Y(NDMAX), FTABLE 
COmON / APPROX / BREAKChriAX), COET ( KTHWAX ) , L, K 
, V(M1AX) 

NDATA ■ 55 
^PITE(20,5) 

5 FCRnAT(' INDY DATAV/' YEAR Y X') 

DO 1 I-L NDATA 

READ(21,4) YEAR, YCI),X(I) 

4 F0R1AT(I4, lX,r7.3, 1X,F2.0) 

WITE(20,2) VEAR, Y(I),X(I) 

2 F0RMAT(I4, 1X,F7.3,1X,F3.0) 

1 CCNTINJE 

C GIVE T>€ ORDER K AND MJT1BER OF INTER'/ALS L 

K ■ 4 

L • 3 

FTABLE -8.00 

C GIVE THE BREAKPOINTS AND CONTINUITY CONSTRAINTS 
BREAK(l) ■ 0, 

BREAK(2) • 7.5 


59 



OF f^OR QOAUTT 


BREflKO) - 33.5 
BR£Ak(4) • 62. 

V(l) «0 
V(2) - 3 
V(3) • 3 

RETURN 

E?® 

SUBROUTI^C M.3<r(BREflK,L,K,V,T,N.KEND) 

CCWUTES TVE KNOT SEQUENCE T AND DI^€NS10N N FROM TVC BREfW=OINT 
C SEQUEhCE BREAK, GIVEN TVE SR_Ih€ ORDER K, THE tSUMBER OF INTER- 
C VALS L. AND TVE NUHBER OF CCNTINUITV CONDITIONS V(I) RT BREAK 

C 

input 

C BREAK Cl) BREAK(L+1) TVE BREAKPOINT ^QUENCE. 

C I TVE NUHBER OF INTERVALS. 

C K TVE ORDER OF TVE SPLItE. 

C V(2), . . . ,V(L) TVE VUtBER OF CONTINUITY CONSTRRINTS RT 

C BREAK(2)....,BRERK(L). 

C 

C)> «***» OUTPUT *> n c*i*m 

C T(l), T(N+K) TVE KNOT SEQUENCE. 

C N....TVE DIMENSION OF TVE SPLINE SPACE OF ORDER K. 

C KEND(1)....TVE INDEX OF TVE LARGEST KNOT EQUAL TO BREAK(I) 

METHOD 

C TVC FIRST K KNOTS ARE SET EQUAL TO BREAKCl). TVE KNOTS ARE 
C TVEN SEQUENCED SO THAT K - V(I) KNOTS ARE AT BREAK(I) WITH 
C KEND(I) EQUAL TO TVE INDEX OF TVE LARGEST KNOT AT BREAK(I). 

C N IS SET EQUAL TO KENDCL) AND TVE LAST K KNOTS T(hH-l), . 

C T(NHO ARE SET EQUAL TO BREAKCL>1). 

C 

INTEGER K,L,N, I. ISTART, ISTOP, KEND(l) 

REAL BREAK(l), T(l) 

C SET TVE FIRST K KNOTS EQUAL TO BREAK(l). 

DO 1 I • 1, K 

1 T(I) • BREAK(l) 

C 

C FIND TVE INDEX KEND(I) OF TVE LARGEST KNOT EQUAL TO BRE»<(I) 
KEND(l) • K 
DO 2 I • 2, L 

2 KEND(I) • KEND(I-l) + K - V<I) 

C 

C SET T(KENDd-l) + 1) -...- T(KEND(D) - BREAK(I). 

DO 10 1 ■ 2, L 

ISTART > KEND(I-l) +1 
ISTOP • KEND(I) 

DO 11 J ■ ISTART, ISTOP 

11 T(J) - BREAK(I) 

10 continue 

N • KEND(U 
C 

C SET TVE LAST K KNOTS EQUAL TO BREAK(L+1). 
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DO 20 I • 1, K 

20 T(N+n • BREflk(L+l) 

RTTIWI 

END 

SUBPOUTirC 0JTNTS(BREPK,V,L,T,N,K,KEND) 

C THIS SUBROUTirC IS TOP OLfTPUTING OLY. IT OUTPUTS PLL 
C CflLLING flRGUrCNTS « F T E R VL2KT HAS BEEN CALLED. 

C 

INPUT AND OUTPUT>ii*<i*** 

C K....TVC SPLirC ORDER 

Cl T>€ MJMBER OT INTERVALS 

C H. ,..T>C DlfCNSION OF TV€ SPLirC SPACE 

C BREAK(l), . . .,BREAK(L+1) TVE BREAKPOINT SEQUENCE 

C V(l),...V(L) T>€ MJriBER OF CONTINUITY CONSTRAINTS AT 

C BREAKC1),...BREAK(L) 

C Ta),...T(N).,..T>€ KNOT SEQUENCE 

C KENDa KENDCL) INDEX OF T>€ LARGEST KNOT EQUAL TO 

C BREAK(1),...,BR£AKCU 

C 

DirCNSION T(i), KENDCl), BREAKCl) 

INTEGER V(l) 

I49ITEC20,40^ K, L, N 

40 rORMPTC//' TVE ORDER K • 13//' TVE ♦ INTERVALS L 13, 

* //' TVE DirCNSION N . 13) 

UPITE(20,41) 

41 FOfWTC//' BREAKPOINTS', T20, ' CONTIhtJirf' CONDITIONS') 
DO 45 J • 1, L 

45 H?ITE(20,42) BREAK(J), V(J) 

42 FORf1ATCFi6.8,T30,I3) 

LPl • L 1 

URITE(20,43) BREAK(LPl) 

43 FCRMAT(F16,8) 

URITE(20,8) 

8 F0RhAT(//' T INDEX') 

hPK > N + K 

icoim . 1 

INDEX • i 

UPITE<20,5) T<1), INDEX 
5 F0RriAT(F16.8, 5X, 13) 

DO 7 J • 2, hPK 

IF (T<J) .EQ. TCJ-D) GO TO 50 

UFITEC20, 12) ICOUNT, KZND(ICOUNT), T(J), J 

12 FORflAT(T30, 'KEND(M3, ')• M3/T1 , F16. 8, 5X, 13) 

GO TO 13 

50 >RITE(20,9) T(J), J 

9 F0RrttT(F16.8,5X,I3) 

GO TO 7 

13 ICOUNT . ICOUNT 
7 CONTINUE 


E3SD 

SUBROUTirC nLAG(IFLAG,N) 
C THIS SUBROUTirC OCCKS FOR 


RETURN 
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C (1) BREAKPOINTS U^lICH ARE NOT STRICTLY INCREASING; 

C (2) TOO rWNY CONTINUITY CONDITIONS; 

C (3) K LARGER THAN 20 

C (4) X VALUES OUT OF RANGE OT 1>C FIRST AND LAST BREAKPOINTS. 

C 

PAf»^fCTER(NhAX-i00, NEm<-200,KTTr^ 

INTEGER V 

COmDN / DATA NDATA, X(NDMAX), V(NDmX), FTABLE 
COTflON / APPROX / BREPK(^r1AX),COE^(K^rH1AX),L,K,V(^^^ 

DO 1 I»1,L 
IPl-I+l 

ir(BR£PK(I) .GE. BREAK(I-kD) GO TO 2 

1 CONTIHJE 

GO TO 4 

2 t^ITE(20,3) BREAKCI),BREAKCIP1) 

3 FCRMAT(/' BREAKPOINTS HLST BE STRICTLY INCREASING.'/ 

» ' BREAKPOINT', ri6. 8, 2X, 'IS NOT LESS THAN BREAKPOINT', 

F16.0) 

IFLAG-1 

4 DO 5 I-1,L 

IF(V(I) .GE. K) GO TO 6 

5 CONTINUE 

GO TO 20 


6 H9ITE(20,7) I,V(n,BREAK(I) 

7 rO^W^TC/' T>€ H-NBER OF CONTINUITY CONDITIONS MJST BE STRICTLY'/ 

» 5X, ' less THAN JVC SPLIhE ORDER K, V( M2, ' ) • ' , 12/ 

» 5X, ' AT BREAKPOINT ',ri6. 8, " IS TOO LARGE, ' ) 

IFlAG-1 

20 IF (K -.GT. 20) GO TO 8 
GO TO 10 
0 i^ITE(20,9) K 

9 FORMATC/' K-',I2, ' IS TOO LARGE.'/' THE ORDER K MUST BE 20 OR', 

♦ ' IXSS. ') 

IFlAG-1 

10 DO 11 I-1,NDATA 

IF(X<I) .L£. BREAKd) .OR. X(I) .GE. BREAK(L+1)) GO TO 12 

11 CONTINUE 

GO TO 14 

12 UPITE(20, 13) I, X(I), BREAKd), BREAkCL'*^!) 

13 FORMATC/' X VALUE OUT OF RANGE. '/' X(',I4, ')-',F16.8, 

* MS NOT IN TVE RANGE BREAKd )- \ FIS. 8/ 

m 5X, 'TO BREAK(LAST)-',F16.8) 
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IFL«-l 

14 IF(N .GT, NDATA) GO TO 16 
GO TO 16 

16 I4?ITE(20,17) N, NDATA 

17 TVC DIhGMSION N-M2, 

* ' IS GREATER THPN TVE SAMPLE SIZEM4, '. ') 

IFLAG-l 

18 KTN-K*M 

IFCNDMftX .LT, NDATA) GO TO 19 
GO TO 22 

19 l^I7i:(20,21) NDMAX,NDATA 

21 F0Rf1AT(/^ OCO PARAtCTER STATEMENT. '/5X, ' NDmK-M5, 

* ' rUST NOT BE LESS THAN TV€ NUMBER OF DATA POINTS', 

» 14,'.') 

IFLAG-1 

22 IF (rm< ,LT. N) GO TO 23 


GO TO 25 

23 U«ITE<20,24) N 

24 FORhATC/' OCCK PARATCTER STATOEKT. '/5X, ' M1AX MJST NOT BE', 

» ' LESS THAN N-',I5) 

IFLAG-1 f 

25 IF(KT>f1AX .LT. KTN) GO TO 26 
GO TO 28 

26 URITE(20,27) KTN 

27 F0PMAT(/' OCCK PARAfETEP STATEMENT. '/5X, ' KTHlAX MUST NOT', 

» ' BE LESS THAN KTN- ' , 14, ' . ’ ) 

IFLAG-1 

28 IFCIFLAG. EO. 0) GO TO 30 
UFITE(20,29) 


* ’ STTPDOMH 0»t0T PROCETD. PROGRflfi ABORTS. ') 

30 RmjRN 

END 
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SUBROOTirC LSTSQl(T,N,k,Q,DIAG,a:OeF) 

QPLLS BSPLVB, B06LV 

c 

comow STflTUtJfT DflTft IS USED. 

c 

C THIS IS A tlODiriCATIOH OF DE BOOR'S SUBROUTI^€ L2RPPR, 

C PRGE 255. IT Ihfaijrs T,N,K, Flf^DS TV€ LEftST SQUARES 
C «=PROXimTIOH TO TVC DRTA USING WOW< RRROyS Q AND DIA6, 

C AND OUTPUTS TVE B-SPLItC COETFICIENTS BCOEF. 

C 

PAfiA»-CTD?(l<m><-20, NDflAX-200) 

REAL BCOEF(N), DIAG(N), Q(K,N), T(l), BIATXCKMPW) 

CCW1CN / DATA / NDATA, X(NDnAX), Y(NDflAX), FTABLE 
C 

DO 7 J-l.N 

BCCCFCJ) • 0. 

DO 7 1*1, K 
7 0(1, J) • 0. 

LEFT • K 
LEFTTK • 0 
DO 20 U*l, NDATA 

C LOCATE LEFT ST X(LL) IN (T(LEFT). T(LEFT+1 ) ) 

10 IF CUJT .EQ. N) GO TO 15 

IF (X(LU .LT. T(LEFT+1)) GO TO 15 

LEFT • L£FT-t-l 
LEFTTK • LEFTTK+1 

GO TO 10 

15 CALL BSPLVB(T,K, 1,X(LL),LEFT,BIATX) 

DO 20 W-1,K 

DU • BIATX(m) 

J • LEFTTKHti 

BCOEF(J) • DU»V(LL) + BCOEF(J) 

I-l 

DO 20 JJ-IT1,K - 

Q(I,J) • BIATX(JJ)»DU + Qa,J) 

20 I • I+l 

CALL ft>FAC(Q.K,N.DIAG) 

CALL BCHSLV(Q.K,N, BCOEF) - , 

RETURN _ , . , 

END 

SUBROUTIhC BSPLVB(T,JHIGH, INDEX, X, LEFT, BIATX) 
calculates the value CF all POSSIBLY NONZERO B-SPLlhCS AT X OF ORDER 

JOLfT-MAX C JHIGH, ( J-f-1) •( INDEX-1 ) ) 

WITH KNOT SEQUENCi: T. 

DE BOOR PAGE 134-135 
PARANETER(JNAX-20) 

INTEGER INDEX, JHIGH, LEFT, I,J,JPl 

REAL BIATX(JHIGH),T(1),X, DELTALt JNAX) , DELTAR(jmX) , SAVED. TERN 
C Dirt?€ION BIATX(JOUT), T(L£FT+JOUT) 

C CURREJTT FORTRAN STANDARD MAKES IT irF>OSSIBLE TO SPECIFY TVC LENGTH 
C OF T AND OF BIATX PRECISELY WITHOUT TVC INTRODUCTION OF OTVERUISE 
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C SLPERTLOUS PDDITICm. PRGLICNTS. 

DPTA J/1/ 

C SAVE J,L£LTAL,DaTAR (VALID IH rOfiTRAN 77) 

C 

GOTO (10,20), INDEX 

10 J-1 

BIATXd)-!. 

ir (J.GE.JWIGH) goto 99 

C 

20 JPl-J+1 

D£LTAR ( J ) -T (L£rT+ J ) -X 
DELTAL ( J ) -X-T ( LEFT+i- J ) 

SPvO-0. 


DO 26 >1,J 

TIRn-BIATXt I )/(D£LTAR( I )+DQ.TflL( JPl-I ) ) 

BIATX ( I ) -SAVEIHDGLTPP( I ) *TOTl 
26 SAVCD-DCLTAL(JP1-I)>»TERM 

BIAT^tJPD-SPVED 
J-JPl 

IF (J.LT.JHIGH) goto 20 

C 

99 RETURN 

SUBROLfTirC BCHFfC (U, NBPNDS, MWW, DIPG) 

C 

C COf^TRUCTS TVC CHCLESKY FACTORIZATION C - L * D • L-TRANSPOSE. 
C SEE DE BOOR P. 256 
C 

I^fTEGER NBPNDS, l^?OU, I, imx, J, JMAX, N 
REPL U(NBPNDS,/«)W), DIAG(^«50W), RATIO 
IF ( r«OW .GT. I ) GO TO 9 

IF (W(l,l) ,GT.0.) 14(1,1) - l./UCl,l) 

Kk. I lRN 

C STORE DIAGONPL OF C IN DIAG. 

9 DO 10 N-l,M90W 
10 DIAG(N) • W(1,N) 

C FACTORIZATION 

DO 20 N-l,h«C*4 

IFCWCl^NJ-^-DIAGCN) .GT. DIAG(N) )GO TO IS 
DO 14 J-1, NBPNDS 

14 W(J,N) - 0. 

GO TO 20 

15 U(1,N) ■ l./W(l,N) 

IMAX • HlN0(NBPNDS^l,hROW - N) 

IF (IMAX .LT. 1) GO TO 20 

JT«< * inAX 
DO 18 I-1,IHAX 

RATIO • W(I+1,N)*W(1,N) 

DO 17 J-l,JMAX 

17 W(J,f^I) • I4(J,N+I) - W(J+I,N)#RAT10 

JhAX - JTIAX • 1 
10 W(I+1,N) ■ RATIO 

20 CONTINUE 
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retlrn 

END 

SUB RC tfTIIC BGHSLV(Wj^NBeWI^Wa^B) 

C SOLVtS T>€ LihCPfi SYSTEM C»X-B OF ORDER NROW FOR X 
C PROVIDED W CONTAINS TVE CHOLESKY FACTORIZATION FOR 7>€ BANDED (SYM- 
C I’CTRIC) POSITIVE DEFINITE MATRIX C AS CONSTRUCTED IN TVE SLfflROUTI^E 
C BOTAC <QUO VIDE). 

C DEBOOR PAGE 258 

INTEGER NBANDS,M?OW, J, JMAX,N,NBND«1 
REAL U(NBANDS,r490U),B(h«?0U) 

IF (h«OW.GT.i) GOTO 21 

B(l)-B(l)*«(l,l) 


C FORWARD SUBSTITUTION. SOLVE L>»Y-B FOR Y, STORE IN B. 

21 NBNDMl-NBANDS-1 

DO 30 N-1,K?0W 

JMAX-MirCCHBNDMl , r«OW-N) 

IF (JMAX.LT.I) GOTO 30 

DO 25 J-1,JMAX 

25 B(J+N)-B(J+N)-W(J+1,N)»B(N) 

30 CONTIME 

C 

C BACKSUBSTITUTION. SOLVE L-TRANSP.X-B>wt(-1)*Y FOR X, STCRE IN B. 

DO 40 N-f'ROU, 1,-1 
B(N)-B(N)*W(1,N) 

JMAX-MIN0( NBNDNl , N50W-N ) 

IF (JMAX.LT.I) GOTO 40 

DO 35 J-1,JMAX 

35 BCN)-B(N)H4(J+1,N)*B(J-»N) 

40 CONTIMJE 

RETURN 

END 

SUBRCUTirE bspLpp (T,BC0EF,N,K,SCRTCH, break, COEF.L) 

C CALLS BSPLVB 

c 

c CONVERTS TVE B-REPRESENTATION T, BCC€F,N, K OF SOTE SPLirC INTO ITS 
C PP-REPRESENTATION BREAK, COEF, L, K. 

C DE BOOR PAGES 140-141 
PARATETER ( KMAX ■ 20 ) 

INTEGER K,L,N, I, J, JP1,KMJ,L£FT,LS0FAR 

REAL BCOEr(N),BREAK(l),COEF(K,l),T(l), SCRTCHCK.K) 

•, BIATX(KMAX),DIFF,FKMJ,SLN 

C DI^ENSION BREAK(L+l),CO£F(k,L).T(fHK) 

LSOFAR-0 

BREAK(1)»T(K) 

DO 50 L£FT*K,N 

C FIND TVE rCXT NONTRIVIAL KNOT INTERVAL. 

IF (T(L£FT+1).EQ.T(LEFT)) GOTO 50 
LSOFAR-LSOFAR+1 
BREAK ( LSOFAR-k 1 )- T ( left-*- 1 ) 

IF (K.GT.l) 

COEF (1 , LSOFAR )• BCOEF ( left ) 


GOTO 9 
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GOTO 50 

STORE TV€ K B-SPLirC COOT'S RELEVWT TO CURRENT KNOT INTERVAL 
C IN SCRTCH(.,1), 

9 DO 10 I-1,K 

10 SCRTCH( 1,1) -BCOErcLErr-K+I ) 

C 

C rOR J-l,./.,k-l, COeWTE TVC K-J B-SPLII^ CCEFT'S RELEViW TO 

C OJRREm KNOT INTERV^ TOR TTC J-TH CCRIVATIVE BY DIOTPENCING 

C THOSE FOR TV€ (J-D5T DERIVATIVE, PND STORE IN 5CRTCH( . , J^l) . 

DO 20 JP1*2,K 
J-JPl-1 
KMJ-K-J 

FKMJ-FLOPT(KMJ) 

DO 20 I-l, KhJ 

DI FT- T ( LXTT+ 1 ) -T ( UETT+ 1 -KHJ ) 

IF (DIFF.GT.0) SCRTCHCI, JPD- 

» ( CSCRTOKI+ 1 , J)-SCRTCH(I, J) )/DIFF)»TKnj 

20 CONTINUE 

C 

C FOR J-0,...,K-1, FIND TVE VALUES AT T(LEFT) OF TVC J+1 
C B- SPLirCS OF ORDER J+1 mOSE SUPPORT CONTAINS TVC CURRENT 
C KNOT INTERVAL FROri THOSE OF ORDER J (IN BIATX), TVCN COMBIfC 
C WITH T»C B-SPLIC COETF'S (IN SCRTCHC . ,K-J) ) FOUND EARLIER 
C TO COfPUTE TIC (K-J-l)ST DERIVATIVE AT TaZFT) OF TV£ GIVEN 
C SPLIIC. 

C NOTE. IF TVC REPEATED CALLS TO BSPLVB ARE TVCUGHT TO GEVCRATE 
C TOO MUCH OVERICAD, TVCN REPLACE TVC FIRST CALL BY 
C BIATXCD-l, 

C and TVC SUBSEQUENT CALL BY TVC STATE^e^fT 
C J-JPl-1 

C FOLLOICD BY A DIRECT COPY OF THE LICS 
C DeLTAR(J)-T(L£FT+J)-X 

C 

C BIATXOl) -SAVED 

C FROM BSPLVB. DELTALCKMAX) AND DELTAR(KMAy) WCXJLD HAVE TO 
C APPEAR IN A DirCNSICN STATETCNT, OF COURSE. 

C 

CAUL BSPLVB(T,l,l,T(U£rT),U£rr,BIATy) 

COGF ( K, LSCFAR) -SCRTCHC 1 , K ) 

DO 30 JP1-2,K 

CALL BSPLVB(T, JPl.2,T(LEm,L£rT,BIAT>:) 

KMJ-K+l-JPl 

SLil-0. 

DO 28 I-1,JP1 

29 SUh-BIATX( I ) «SCRTCH( I , KMJ ) +SUT1 

30 C0EF(KMJ,LS0FAR)-5UM 

50 CONTINJE 

L-LSOFAR 

RETIFN 

END 

SUBROUTIC ERRL2KFTAU, ERROR, ERRDF, SSE, MSE) 
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CflU-S SUBPROGRPn PPVfiLU( INTI3?V) 

THIS SUBROUTlrC CClWrZS THE ERROR SS AND MS. IT IS A 
MODIFIED VERSION OF D£ BOOR'S SUBRO/TlrC L2ERR, PAGE 261. 
MSE IS T>€ OtnrVTED fEAN SQUARED ERROR. 

PARA^'tTEP<l♦1AX-100.^®t1AX-200,t<T^f1AX-2000) 

INTEGER ERRET, V 

f?EAL FTAU(l), ERRCR(l), MSE. Y, X, BREAK. COEF 
DirEHSION FTAUCNDATA), ERROR(NDATA) 

COmX / DATA / NDATA, X(NDMAX). Y(NDMAX). FTABLE 
COmOH / APPROX / BREAK(MIAX), COEFCKTrflAX) . L. K 
* , V(ftlAX) 


SSE-0. 

DO 10 LL-1, NDATA 

FTAU(LL) • PPVALU(BREAK,CC€F.L,K.XajL).0) 

ERROR(U-) • Y(LU - FTAU(U-) 

10 SSE • SSE + ERR0R(U.)*»<2 
MSE • SSE/ERRDF 

return 
END 

REAL FUNCTION PPVALUCBREAK.COCF.L.K.X. JDERIV) 

CALLS 'INTERV' 

CALCULATES VALUE AT X OF JDERIV-TH DERIVATIVE OF PP FCT FROM PP-TO=R 
INTEGER JDERIV.K.L, I.M.NDLiTTY 
REAL BREAK(L).COeF(K,L),X, FTtUDR.H 
PPVALU-0. 

FmjDR-K-JDERIV 

DERIVATIVES OF ORDER K OR HIG«R AREARE IDENTICALLY ZERO. 
IF (FrWDR.LE.0) GOTO 99 

FIND into: I OF largest BREAKPOINT TO TVC LETT OF X. 
call INTERV<BREAK.L,X. I.NDOTfY), 

EVALUATE JDERIV-TH DERIVATIVE OF I-TH POLYNOMIAL PIECE AT X. 
H-X-BREAKa.'i 
DO 10 M-K. JDERiV-fl.-l 

PPVALU* ( PPVALU/FmJDR ) *+++CO£F ( M, 1 J 
10 FmjDR-rmjDR-i 

99 return 

END 

SUBROUTI^C INTERV(XT,LXT,X.L£rT,rffTj«) 

C C0rt=UTES L£FT-MAX(I,1.L£.I.LE.LXT.AND.XT(I).L£.X) 

C DE BOOR PAGE 92 

INTEGER LEFT, LXT.J-TLAG. IHI, ILO, ISTEP, MIDDLE 
REAL X.XT(LXT) 

DATA ILO/1/ 

C SAVE ILO (A VALID FORTRAN STATEMENT IN 1>€ ^€W 1977 STANDARD) 

IHI-ILO+1 
IF (IHI.LT.LXT) 

IF (X.GE.XTCLXT)) 

IF (LXT.LE.1) 


GOTO 20 
GOTO 110 
GOTO 90 
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c 

20 

C 

C 

31 


35 


C 

40 

41 


45 

C 

C 

50 

C 


GOTO 40 
GOTO 100 

wswNOW X.LT.XT(ILO). DECREASE ILO TO CAPTURE X. 


ILO-LXT-1 

IHI-LXT 

IF (X.GE.XT(IHD) 
IF (X.GE.XT(ILO)) 


IS7EP*1 
IHI-ILO 
ILO-IHI-ISTEP 
IF (ILO.LE.l) 

IF (X.GE.XT(ILO)) 
ISTEP-ISTEP^ 

ILO-1 

IF (X.LT.XT(l)) 

X.GE.XTT(IHI). 

ISTiP-1 
ILO- INI 
IHI-ILOISTEP 
IF CIHI.GE,LXT) 

IF CX.LT.XT(IHI)) 
I5TIP-ISTEP*2 

IF (X.GE.XTOXn) 

IHI-LXT 


GOTO 35 
GOTO 50 

GOTO 31 

GOTO 90 
GOTO 50 

INCREASE IHI TO CAPTURE X. 


GOTO 45 
GOTO 50 

GOTO 41 
GOTO 110 


XT(IL0).L£.X.LT.XT(IHI). narrow n-C INTERVAL- 
MIDDLE- (ILCKlHI)/2 

IF (MIDDLE. EQ. ILO) GOTO 100 

NOTE. IT IS ASSLM3) THAT MIDDLE- ILO IN CASE IHI-ILO+1. 

IF (X.LT.XTCfllDDLE)) GOTO 53 

ILO-MIDDLE 

GOTO 50 


53 IHI -MIDDLE 

GOTO 50 


C<e#»SET OOTPUT and RETURN. 


90 

rTLAG— 1 
LEFT-1 




RETURN 

100 

fTLAG-0 

LEFT-ILO 

Rrn«H 

110 

hFLAG-l 

lett-lxt 

RETURN 


END 

SUBF?OLmr€ CNTRSTCI, JDERIV, N, K, L, T, BREAK, KEND, 

* BRT, BLE, DCOEF, CTRAST) 

CALLS BCONT 
C 

C finds T>C CCrfTRPST CCETFICIEKTS FOR TESTING CONTINJITV OF THE 
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C JICRIV-TVI ECRIVATIV^ OF THE SPLlrC njNCTIOM AT BREPK(I). 

C 

C <nm > n < nr» INPUT <DMo«ot 

C I NLflBER or INTEPVPLS 

C T(1),...,T(N+K) T>C KNOT SEQUE3<£ 

C I TVE INDEX or TV€ BREAKPOINT OT INTEREST 

C BREAK(1),.,.BREAK(L-M) T>€ BREAKPOINT SEQUENCE 

C JTERIV NOTfCGATlVE INTEGER GIVING T>€ ORDER OF TVE DERI- 

C VATIVE TO BE EVALUATED 

C KENDU), . . .KZND(L) INDEX OF THE LARGEST KNOT EQUAL TO 

C BREAKd BREAKCU 

C N DirENSION OF SPLI^€ SPACE 

c K.... ORDER or SR.irc 

C BRT, BLF, KOEF...WORK ARRAYS OF LENGTH N 
C 

C<uMnMO» OUTPUT 3WOWDW 

C CTRAST(1),...,CTRAST(N) TVE CONTRAST CCEFTICIENTS USED TO 

C TEST CCrrriNUITY OF TVE JDERIV-TH 

C DERIVATIVE AT BREAK<D 

C 

Qjiofomtotoi n E T H 0 D 

C TVE FUNCTION SUBPROGRAM BCOMT IS USED TO COTFUTE THE VALUE OF 
C THE LEFT AND RIGHT LIMITS OF TVE JDERIV-TH DERIVATIVE OF 
C RELTVAffT B-5R.INES AT BREAK(I). 

C 

INTEGER KEND(l) 

REAL BRTCl), BLF(l), CTRAST(l), T(l), 

» BREAKCl), Dcoerci) 

DO 20 JJ • 1, N 
20 XOEF(JJ) • 0. 

DO 10 J • 1, N 
DCCEFCJ) • 1. 

COmjTE VALUE FOR RIGHT CONTINUITY 

IF (MEND(I)-K+1 ,L£. J .AND. J .L£.K£NDa)) GO TO 30 
BRT(J) ■ 0. 

GO TO 40 

30 BRT(J1 • BCONT(T,DCOEr,N,K,BREAKCI),KEND(I), 

JDERIV) 


CCWJTE VALUE FOR LETT CONTINUITY 

40 IF(KEND(I-1)-K+1 .LE. J -AND. J .L£. KENDd-t)) GO TO 50 
BLF(J) • 0. 


GO TO 60 

50 BLFCJ) • BCONT(T,DCOer,N,K,BREAK(I), 

♦ KEND(I-l), JDERIV) 

C 

CCmjTE DIF7ERENCE CF TVE LETT AND RIGHT VALUES 
60 CTRAST(J) . BRT(J) - BLFCJ) 

X0E7(J) • 0. 

10 CONTINJE 
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END 


RETURN 


REAL RtCTION BCONT(T,BCCeF,N,K,X, I, JDERIV) 

CALCULATES VALUE AT X OF JDCRIV-TH DERIVATIVE OF SPLIhC FROM 
THIS IS A MODIFIED VERSION OF DE BOOR'S SUBROUTirC BVALUE, 
PAGE 144. THE 0K.Y DIFFERENCE IS THAT T>€ LEFT-HAND KNOT 
INDEX I IS I^«JTED RATVCR THAN FOUND IN INTERV. CONSE- 
QUENTLY, Lire 10 IS MODIFIED TO IM=>UT I AND LI^eS 710 AND 
720 ARE OMITTED. T>€ PURPOSE IS TO ALLOW EVALUATION AT 
B?EAKPOINTS WITH LEFT (0R RIGHT) CONTINUITY. 
PARAFETER(KMAX-20) 


B-REP. 


INTEGER JDERIV, K,N, I, ILO, INK, J, JC, JCMIN, JCMAX, JJ,KJU,KMl,rTLAG 
• ,rMi 

REAL BC0£F(1),T(1),X, AJ(KNAX),DL(KMAiX),DR(KMAX),FKMJ 

DIhENSION T(NHO 

BCONT-0. 

IF (JDERIV. GE.K) GOTO 99 


ww IF K-1 (AND JDERIV-0), BCONT-BCOEF( I ) . 
KMl-K-l 

IF (KJ11.GT.0) GOTO 1 

BCONT-BCOEF(I) 

GOTO 99 


STORE TVE K B-SPLIhC COEFFICIENTS RELEVENT FOR TVE KNOT INTERVAL 
(T(I),T(I-t-l)) IN AJ(l), . . . ,AJ(K) AND COtWTE DL( J) •X-f( I+l-J) , 
DR(J)-T(I+J)-X, J-1,..,,K-1. SET ANY OF THE AJ NOT OBTAINABLE 
FROM DWr TO ZERO. SET ANY T.S NOT OBTAINABLE EQUAL TO T(l) OR 
TO T(NH<) APPROPRIATELY. 

1 JCMIN-1 
IfK-I-K 

IF (IN<.GE.0) GOTO 9 

JCMIN-l-irK 

DO 5 J*1,I 

5 DL(J)»X-T(I+1-J) 

DO 6 J*I,KM1 

AJ(K-J)-0. 

6 DL(J)-DLd) 

GOTO 10 

8 DO 9 J-1,KM1 

9 DL(J)-X-T(I+1-J) 

10 JCMAX-K 
N1I-N-I 

IF (M1I.GE.0) ■ GOTO 18 

J01AX=K+N1I 

DO 15 J-1,J01AX 

15 DR(J)-T(I+J)-X 
DO 16 J-JCMAX,KM1 

Aj(J+l)-0. 

16 Dfi(J)-DR(JOlAX) 

GOTO 20 

18 DO 19 J-l.KJIl 
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19 DP(J)-T(I+J)-X 
C 

20 DO 21 JC-JCHIN,JCmK 

21 AJ(JC)«BCOerciM<tJC) 

c 

C «» DirFERENCE T>€ COEFriCIENTS JDERIV TlrCS. 

IF (JDGRIV.EO.0) GOTO 30 

DO 23 J-1,JDCRIV 
KHT-K-J 

FknJ-FLOAT(Km) 

2 L 0 -KnJ 

DO 23 JJ-l,Kjnj 

ftJtJJ)-{(AJ(JJ+l)-AJ(JJ))/(DL(ILO)+DR(JJ)))*TKnJ 
23 ILO*ILCHl 

C 

C COrfVTE VI^CE X IN ( T( I ) , T( I+l ) ) OF JDERIV-TH DERIVATIVE, 

C GIVEN ITS RELEVENT D-SPLirC COEFS IN AJ( 1 ) , . _ , AJ(k-JDERIV) . 

30 IF (JlCRIV.EO.Km) GOTO 39 

DO 33 J-JI€PIV+l,Km 
kWJ-K-J 
ILO-WU 

DO 33 JJ«l,Knj 

AJ(JJ)-(AJ(JJ+1)»DL(IL0)+AJ( JJ)*DR(JJ))/(DL(ILO)+DR( JJ) ) 


33 

ILO-ILO: 


39 

BCONT-AJ(l) 


99 

END 

RETLRN 


C THIS IS FOR 1 DF HYPOTVCSES. 

SUBROUTI^€ SSHYP(BCOEF,CTRPST,W,NBM)S,N,PVAR,VRR, 

* SSH, rOH, HDF) 

CPLLS FORSUB 
C 

C FINDS TV€ VARIRhCE OF CONTRAST AND TVC MS FOR TESTING THAT 
C T>€ CONTRAST IS ZERO, 

C 

C LINV. , _TV€ INVERSE OF L OBTAirCD FROM B C H I N V 
C CTRAST, . . . TV£ CONTRAST VECTOR OBTAirCD FROM C N T R S T 
C BCOEF TVE S-SFLirC COEFFICIENTS 

C W. . . .TVC MATRIX FROM B C H F A C CONTAINING D-INVERSE 
C NB^^^^DS. , , E(XIALS X 

C N Tl^'NjriBER OF ELDCNTS IN T>€ CONTRAST VECTOR— 


C ALSO TV€ DI^f>eION OF T>t SPLI^€ SPACE 

C PVAR W0»< VECTOR OF LENGTH N EQUAL TO T>€ PRODUCT 

C W(1,.)*A, I.E. D-INSNtLiNVNCTRAST 

C 

C#*#***OUTPUT****#» 

C VAR COEFFICIENT OF SIGMA-50UARED IN T>C VARIANCE OF T>€ 

C CONTRAST. I.E. 1>€ PRODUCT CTRAST-TRANSPOSE*LlNV- 

C T7V^SPbSE*D-INV^INV<:TRAST 

C SSH, MSH. HDF TVC SS, MS, AND DF FOR TVC HYPOT)-ESIS 

C 
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C»^****mCTH 0 D*‘^*** + 

C Th€ PPODUCT LINV<TRAST IS OBTAirCD T^€N PRO'ULTIPLICD 
C D-INV T>€N THPT RESULT IS PREMJLTIPLIED BY ( LINV^CTRAST )- 
C TRPHSPOSE 
C 

IhfTEGER KDT 
REAL NLM, M9H 

REAL CTRAST(l), W(NBAHDS,N), P>/AR(N) 

REAL BCOEFd) 

Nun- 0. 

DO 3 

3 NLTI • NLM + (CTRASK 1 1) -tBCOGF< 1 1) ) 

CALL FORSUB(U,CTRAST,NB(^DS,N) 

DO 1 II-1,N 

1 PVAR(II) ■ W(1,II)KTRAST(II) 

VAR • 0, 

DO 2 II-UN 

2 VAR • VAR + CTRAST(II)*FVARai) 

SSH ■ <NLm<«2)/VAR 

MSH - SSH 
HDT • 1 

RETURN 

END - 

SUBROUTI^E ^0RSUBCU,AA,NBA^4DS,^«)U^ 

SOLVES LV-AA FOR Y AND STORES IN AA 

jif*m***IfHPUT<'***** 

U. . . A MATRIX FED IN FROM B C H F A C AND CONTAINING IN ITS ROWS 
THE DIAGONALS OF A P, D. SYMhCTRIC MATRIX C 
NBANDS. . .!>€ BANDWIDTH OF C 
NROW. . ,Th€ ORD or C 

AA. . .TVE VECTOR OF LENGTH l^?OW CONTAINING T>€ RIGHT HAND SIDE 
»*#*** 0 UTPUT**<f**^ 

AA...TVC VECTOR OF LENGTH I^POW CONTAINING THE SOLUTION 
*j#***#nETHOD^*^**>*f 

THE FORWARD SUBSTITUTION RaJTI^€ FROM. DGBOOR'S BCHSLV IS USED 

REAL WCNBANDS, hFOW), AA(hFOW) 

IF CNROW.GT, n GO TO 21 ■ 

AAC1)-AA(1)*W(1, 1) 

RETURN 

21 NBNDm=NBANDS-l 

DO 30 N«l,NPOW 

JMAX-MIN0( NBNDMl , r4?OW-N ) 

IF (JMAX.LT. 1) GO TO 30 
DO 25 J-1,JMAX 

25 AA(J+N)«AA(J^)-W( J-d,N)<^AA(N) 

30 CONTINLE 
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EDSD 

SUBROimrC REXhCT(KDiD,W<)T,N,K,L,T,V, break) 

C RELABELS tVC KNOT SEQUENCE T( ) BY OMITTING 1>C LEAST 
SIGNIFICANT KNOT, BREAK(KNOT) 

KEND(I). . ,T>€ INDEX OF T>C LARGEST KNOT EQUAL TO BREAKd) 
KNOT, . . INDEX OT T>€ BREAKPOINT TO BE OMITTED 
N. . .DIl^DiSlON or TVC (CLD) SPLI^€ SPACE 
K. . .ORDER OF T>€ SPLIhC 
T... KNOT SEQUENCE 

V(I).. NUMBER OF CONTINUITY CONDITIONS AT BREAK(I) 

BREAK. . .BREAKPOINT SEQUENCE 


*3Mi<nMoMo»D»oMD» C LrrPUT < C3Mc »JM P<Q^^ 

N. .. DIMENSION OF (h€W) SPLIhC SPACE UITH BREAK(KNOT) OMITTED 
T(l>. ,T(N), . . (hEW) KNOT SEQUENCE WITH BREAK(KNOT) OMITTED 


SINCE BREAKCKNOT)-T(KEND(KNOT-l)>l)-. _ -TCKEND(KNOT) ), ut 
RELABEL ALL T^S BEYOND. 


DI/-ENSION KEND(l), T(l), BREAK(l) 
INTEGER V(l) 

I1-KEND(KN0T-1)+1 

I2«KEND(KN0T)'H 

jl»N+K-I2+l 
DO 1 KT-1,J1 
Kl-KT-1 

1 T(I1+K1)-T(I2-H<1) 
N-N-CK-V(KNOT)) 

DO 2 II*KNOT,L 
BREAKai)=BREAK(II + l) 

IF(II .EQ. U GO TO 2 

V(II)-V(II+l) 

KtND<II)-KEND(II-l )+K-V(II) 

2 CONTiriJE 

L-L-1 


RETURN 

END 

SUBROUTirC STDeP9(W,BC0er,K,N,L,MSE,BB,AA,SE,LlNV) 
BCHINV AND MATVEC 


THIS SUBROL^TI^€ COn=UTES T>C STANDARD ERRORS CF THE B-SPLIhE 
COEFTICIENTS AND OUTPUTS THEM, 


REAL W(k,N),BCCCF(N),«SE, BB(N,N) , SE(N) , LINVCN, N), AA(N, N) 
CALL BCHINV(W,K,N,LINV) 

WRITEC20, 10) L,K 

10 FORMATC///' PROCEDURE TERMINATES WITH L-M3, ' AND K-M3// 
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* ' N CC€F S.E.') 

DO 11 II-1,N 
DO 11 JJ-l,N 

11 JJ) 

CALL MATVEC(N,N,N,BB,Lirfv',AA> 

DO 13 II«1,N 

SE ( 1 1) • SQR T ( AA (1 1 , 1 1 ) *liSE > 

UI?ITr(20,12) II,BCOEF(II),SE(II) 

12 rORMAT(I3,2F16.0) 

13 CONTINUE 

RETIFN 

END 

SUBROUTirC BOHIW CW, NBANDS, r«DW, INV) 

C FINDS L-INVEPSE UfCRE L IS T>€ LOWEP TRIANGULAR NATRIX 
C IN T>C OOLESKY FACT0RI2ATICN OF THE BANDED SVtfCTRIC P.D. 

C MATRIX C AS CONSTRUCTED IN T>€ SUBROUTirC B C H F A C. 

C SEE DE BOOR, P. 256 
C 

C«3ia»aK« INPUT 

C hROW.....IS TIC ORDER OF T>E MATRIX C. 

C NBANDS IS nc BANDWIDTH OF C. 

C W. ....CONTAINS THE CHOLESKY FACTORIZATION OF C AS OUTPUT 
C FROM SUBROUTIIC B C H F A C WITH ROWS 2 THROUGH NBANDS--1 
C CONTAINING TC NOTf-ZERO AND NONHJNIT DIAGONAL ENTRIES 

C OF L, 

C 

C ^MoM otoo OUTPUT 

C INV TVC INVERSE OF L. 

C 

CjWuMdW method 

C TC LIICAR SYSTEM L*L-INVEPSE • IDENTITY IS SOLVED FOR 
C L*INVERSE BY SUCCESSIVELY FINDING THE COLUTTC OF L-INVERSE 
C USING T>€ FORWARD SUBSTITUTION ROUTlrC IN B C H S L V. 

C 

INTEGER NBANDS, N?OW, J, JMAX,N, NBNDMl 
REAL W(NBANDS,M50W), INVChROW, rPOW) 

IF ( hROW .GT. 1 ) GO TO 21 

INVa,l) ■ 1. 

RETLf?N 

C 

C STORE TVC IDENTITY MATRIX IN INV 
21 DO 10 J-l,rPOW 

DO 10 i«i,r«ow 

IF (I ,EQ. J) GO TO 20 

INV(I,J) • 0. 

GO TO 10 

20 INV<I,J) . 1. 

10 CONTINUE 

C 

C NOW USE FORWARD SUBSTITUTION FROM 8 C H S L V. 
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NBNDm • NBflNDS - 1 
DO 40 J-l,r^!OW 
DO 30 ,Nrl,^f?OW 

■ MIN0(rHBNEni,lf!pU-«) 

ir (JMflX .LT. 1) (50 TO 30 

DO 25 I • 1,JMRX 

g 1W(14N,J) . INV(I-f«.J) - W(I+1,N)*INV<N,J) 

30 CONTIMJE 
40 CONTIHJE 


RETURN 

SUBROl/TIhC MRTVEC(N,M1,M,X,Y,Z) 

C riNDS T>€ WTRIX OR VECTOR Z l^ICH IS 1>€ PRODUCT 
C X*Y 14HERE X IS PHD Y IS r««wV 

C ' ' 


3 

2 


REPi. X(N,M1), Y(hn,ri), 2(N,M) 
DO 1 r - l’,N 
DO 2 J • l.H 
Z(I,J) • 0. 

DO 3 


Z(I,J) - 2(1, J) + X(I,K)»Y(K, J) 
COmiHJE 


1 CONTINUE 


END 


RETURN 


♦ 
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ORDER REDUCTION PROGRAM LISTING 


PROGRfifI >«=LT 2 (ir^\n',OLrrPl/T,TflPE 6 -OUTPlJT.TPPE 20 ,TflPE 21 ) 
C STEPDOm rop reducing spline order (FOR ftj_ INTERVft_S 
C SIMULTANEOUSLY) tWiLE KEEPING T>€ KNOTS FIXED AND ASSUlING 
K-2 COKTINUITV CONDITIOT€ 


NDflAX IS AT least TVC SAff=LE SIZE, NDATA. 
rtlAX IS AT LEAST N. MITH MAXiriil CONTINUITY CONSTRAINTS, 
N-L+K-1. WITH NO CONTINJITY CONSTRAINTS, N-L-K. 
KTM1AX IS AT LEAST K*N. 

PARAfCTER (rm<-100,r€t1AX-200,KTTriA><-2000) 
real BC0EF(^m<),O(KT^m<).DIAG(KT^f1AX),T(ND^1AX) 

» ,LINV(KTTm<),DCOEr(r^ 1 AX),BRT(hriAX),BLF(rf 1 AX) 

» ,AA(ht1AX),VAR(NlAX),B(rt1AX),C(N1AX),ATRP(rriAX) 

» ,F<^E)mX),ERRC)R(NE»1AX),MSH,MSE,SE(^r1AX) 

• ,KMAT(ft1AX,r#lAX),FB<hm<) 

* ,lAAR(l^1AX),CC(ft1AX),CT(rt1AX),CTRAST(I<1AX) 

INTEGER ERRDF,HDF,V,KEND(hMAX) 

COtnON /DATA/ NDATA,X(NDriAX),Y(NDnAX),FTABLE 
COrtlON /APPROX/ BREAK(r#1AX).COeF(KTTm<).L,K,V(M1AX) 

ICOLNT-0 
C ENTER DATA 

CALL DATiaCOUKT) 

C GET TVC KNOT SEQUENCE 

CALL VL2NT(BREAK,L,K,V,T,N.KEND) 

C preliminary OUTPUT 

CALL OUTNTSC BREAK, V,L,T,N,K,KEND) 

OeCK irCUT DATA 
IFLAG-0 

CALL FLAG(IFLAG,N) 
irciFLAG .EQ. 1) GO TO 25 

CALL PSEUDO 

IEND-0 

LMl-L-1 


C WE WILL TEST THAT TVC K-l-ST DERIVATIVE IS ZERO IN ftj. INTERVALS 

C GET TVC LEAST S(3uARES FIT, I.E., TVC B-SPLIVC COETFICIENTS. 

1 CALL LSTSQ1(T,N,K,Q,DIAG.BC0EF) 

C LSTSQI CALLS BSFLVB, BOCAC, AND BCHSLV. 
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C GET SSE AND fiSE 

EPPOr-MDAT^ 

CPLL B5PLPP(T,BC0Cr,N,K,DIPG,BREPK,C0Er,L) 

CPLL EPRL21<r,O?R0R,ERRDr,SSE,MSE) 

C O^ra-21 CALLS PPVALU WICH CALLS INTERV, 

IFCIEND .EQ. 2) GO TO 30 

LPl-L+l 
DO 10 I-1,LP1 

10 rBd ) -PFM: jLUC break, COET, L, K, BREPK( I ) , 0 ) 

C)Mo»>r»Tfc**Tr»<nfc R.0T5 

CALL Ihf"OPLT(0,NDATA,K,l,r,l,0. ,62,,74. , 150. , 1. , 

* 9,9HINDY DATA, l,lhnr,0,5.,4. , .75, .75) 

CALL IrrOPLT(0,NDATA,X,l,Y,l,0. ,62. ,74. , 158. , 1. , 

* 9, SHINDY DATA, 1, IKY, 22, 5., 4., .75, .75) 

CALL irrOPLT(l,Lm,BREAK(2), l,rB(2), 1,0. ,62.,74. , 158. , 1 . , 

* 9, SHINDY DATA, 1,1HY,1,5.,4.,.75, .75) 


KNl-K-1 

W^-K-2 

Kra-K-3 

irclEND.EQ.l) GO TO 8 
ir(V(2).EQ.KM2) GO TO 12 
4?ITE(20,15) 

* //' T>C SMOOTHEST SPLirC OF ORDER K-M2,2X, 

* ^WITH MAXIMUM CONTIMJITY CM2,2X^, 

» 'HAS SSE-\ri6.8,2X, 'AND MSE-',ri6.8) 

IF(K.EQ.l) GO TO 8 
i^ITE(20,l6) K,Krt3,Km,KM3 

16 F0RMAT(x^ CAN ORDER K- M2, 2X, 'WITH SUB-MAXIMLM CONTiriJiTY C 

* I2,2X^, 'BE REDUCED TO ORDER K-M2,2X, 

* 'WITH MAXlriJM CONTiNUtlY C', 12, ' ’') 

DO 19 II-2,L 

10 V(II)-K-2 

call VL2NT(BREAK,L,K, V,T,N,K£M)) 

GO TO 1 


C TEST FOR L0U£R ORDER WITH TV€ HYPOTHESIS MATRIX KMAT. 

12 DO 80 JJ-1,N 

80 XOEFCJJ) - 0. 

DO 2 III-1,N 

XOCTdll) • 1. 

DO 81 II-1,L 

KliATdl, III)-XONT(T,XOEr,N,K,BREAKdI),KEND(II),Km) 
M-dI-l)»N+III 

81 CT(M)-KWT(II, III) 

XOEFdll) -0. 
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2 COfTlNJE 

CflLL SSHYP2(BCC€r,a,0,K,L.N,(¥l,DIAG,Vf«,SSW,l1SH,Hr)F,B,C 

C SSHYP2 CPLLS FORSUB W ri^TVEC. 

TOPTIO-liSH/tlSE 

IF(FTOTIO.G£:.rrABL£J GO TO 5 
UFITE(20,3) 

3 FORMPK//" YES. M 

UFITE(20,31) K,Kro,rrftEL£,FRPTIO,SSE,riSE 

31 FORflPK' FOR K- M2,2X, 'ATSD CM2,2X/5X, 

# 'FTRBLE VRLUE • ' , Fl6. 0, SX, 'OBSEFMD F-\Fl6.0 

It /5X, 'SSE-',F16.0,2X, 'MSE- \ F16. 0) 

CPLL V12rFr(BREPK,L,K,V,T,N,KEND) 

GO TO 1 

5 IEND-1 
UPITE(20,6) 

6 FORMPTC//' NO. ') 

145ITE(20,31) K,KK3,rTABL£,rRRTIO,SSE,f1SE 
14?ITE(20,32) 

32 FORMPTC//' PROCEDURE TERMINATES. 

DO 71 II-2,L 

71 

CPLL VL2NT<BREPK,L.K,V,T,N.KZND) 

GO TO 1 

C PRINT RESULTING CCCJTICIENTS AND STANDARD ERRORS. 

0 URITE(20,13) L,K,kW 

13 FORMAT (///' PROCEIURE TERMINATES WITH L- M2, ' ; K-M2, CM2 

* //' N COCT ST. ERR.') 

CALL STDERR(Q,BCOCr,K,N,L,MSE;DIAG,AA,SE,LINV) 

C STDERR CALLS BCHINV AND MATVEC. 

IF(K .EQ. 1) GO TO 25 
IEND-2 

K-K-l 

DO 40 II-2,L 
40 V(II)-K-1 

CALL V12NT(BREAK,L,K, V,T,N,KEND) 

GO TO 1 

30 H?ITE(20,29) KM1,KM3,SSE,MSE 

29 FORMAT <//' <u»*3»jMi****«54i*jMt** FLRTl-CR IhFORMATION / 

m 5X, TOR K-M2,2X, 'AND CM2,2X/ 

« 5X, 'SSE*',F16.0,2X, 'MSE»',F16.0) 

CALL CALPLT(0.,0.,999) 

25 STOP 


END 
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SUBPOOTIhC SSHVP2(BCCCr,CT,W,NBPl^,f<0h,N,fl 
* ,PVW,VPR,S9H,f^,HDT,B,C,f=iTRP, WV<«,CC,CTR^ 

INTEGER Hi:r,NCCM,N,NBfiND5 
REA- SSH 

REA. PVPR(rCOH,N),A(NCON,N),W(NBPNDS,N),CT(l) . 

REA VflPcr<OrS,NCON), B(N), C(NCON), f^TRP(N,NCON) 

REA WVPR(NCON,MCOh), CTWSKNCOrH, N) , BCOETCN) , CC(N) 
DO 1 I*i,NCOS 
DO 2 JJ-1,N 

CCaj)-CT(M) 

2 CTRASTd, JJXCCJJ) 

CAL rORSUB(W,CC,NBPHr)S,N) 

DO 3 J-1,N 

3 A(I,J)-CC(J) 

1 COHTIMJE 

DO 4 II-1,NCCM 
DO 4 JJ-1,N 

4 PVARdI, JJ)-WCl, JJ)^(II, JJ] 

DO 5 I-i,N 

DO 5 J«1,NC0N 

5 ATRPd, J)-A<Jd) 

CAJL ^PT\^C(f<ON,N,^<C^I,PVflR,RTT^,Vf^^ 

DO 6 I-l.NCON 

MI-rCCN-I^l 
DO 6 J-l,m 

6 WVRRd, J) 

CAL 0>f"A (WVRR,NC0N,NC0N,DIA) 

CAL flPTVEC(NCON,Nd,CTPRSf,BCOEr,B) 

CAJ. ^oPSUB<t.^w,B,^co^l,r<:or^) 

DO 8 Jd^NCON 
3 C(J)»WVPR(1, J)*B( J) 

CAj. nATVEC(i,Ncas,i,a,c,ssH) 

f1^«SSH/NC0N 

HDT-NCON 

RETLFN 

END 
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